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PREFACE 


This is Volume II of a three-part report which describes a computer program 
for the prediction of Solid Propellant Rocket Motor Performance. The computer 
program described herein will be refered to as the SPP program. 

Volume I of this report describes the engineering analysis which was used 
in developing this computer program. 

Volume II of this report is a programming document of the computer program 
which was developed under this contract. It includes a subroutine-by subroutine 
description of all o:c the elements of the SPP program. 

Volume III of this report is a Program User's Manual which describes the 
input necessary to execute the SPP computer program and the information required 
to interpret the output. A sample case is also included. 
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1.0 Introduction 

This document, Volume II, provides a detailed description of all of the 
elements used in the Solid Propellant Rocket Motor Performance Prediction (SPP) 
computer program. The goal in developing the SPP code was to develop a computer 
ized analysis technique that could predict solL’ propellant rocket motor delivered 
specific impulse to within +2% and delivered thrust and total impulse to within 

ts%. 

The contents of this volunu are not intended to supply the reader with a 
complete description of the SPP computer code. Instead, it is intended to give, 
along with the information supplied in Volumes I and III, the necessary information 
about the internal functions of the SPP code to allow the educated user to under- 
stand and/or mod if v the basic internal workings of this computer program. 

The SPP code consists of five basic computational modules plus a master 
control module which selects via user input which of the computational modules 
are to be executed. The approach used in developing this code has been to 
supply the user with the maximum amount of flexibility consistent with ease of 
use. 

The five basic computational modules used in the SPP code are described 
very briefly below it. Table 1-1. A more detailed description of these modules is 
presented in Volume T and to a lesser extent in Volume III. 


Module 

Description 

0DE 

Calculates the reference or Realized performance 
of the solid propellant rocket motor. 

BAL 

Calculates the grain shape and internal ballistics 
parameters as a function of time. 

0DK 

Calculates the nozzle performance considering the 
degradation in I due to finite rate chemical kinetics 

Sp 

TD2P 

Calculates the nozzle performance considering the 
degradation in performance due to two dimensional 
and two phase flow effects. 

TBL 

Calculates the nozzle performance considering the 
degradation in performance due to a viscous turbulent 
boundary layer. 
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Volume I consists of the engineering analysis incorporated in the SPP code 
while Volume III describes the input necessary to execute the SPP code. The 
user's manual (Volume III) and this volume are considered to be the prerequisite 
amount of informat.on needed to successfully understand and/or modify the SPP 
program. 

Section 2 of this volume describes the overall flow of the SPP program along 
with the data communicated between the various computational modules. The over- 
lay structure and labeled common blocks used by each subprogram are detailed in 
section 3 while the linkages and files used by the SPP code are described in sec- 
tion 4. A complete subroutine by subroutine description of this code is included 
in section 5. Section 5 is orgainzed by overlay structure and includes a description 
of the main routine of each overlay followed by a description of each subprogram 
in that overlay in alphabetical order. 

Liberal use of existing documentation, References 1 through 6, about the 
'ndividual elements of the SPP code has been used in this volume. The interes- 
ted reader would be well advised to obtain these documents if he does not already 
have access to them. 




2,0 SPP Program Flow Charts 

This section cf Volume II, the Computer Program Description, is intended 
to give the user of the SPP code a general description of the overall flow of the 
computer program. Also described here is how data is transmitted between modules 
when the many options which are available in the SPP code are executed. 

Figure 2-1 shows the overall logical flow of the SPP code. The solid lines 
in this figure show the mandatory flow while the dashed lines indicate the optional 
paths the computer program can take. In the case that more than one optional 
path can be taken, a diamond sign,0 , is used to indicate the alternate paths. 

The boxes on the right hand side of Figure 2-1 show the action taken by 
the SPP code in response to Input directive cards and Iso to the selection of 
loss modules to be executed. The items on the left band side of this figure show 
the alternate logical execution paths which are determined by user input. 

Figure 2-2 shows the actual execution sequence taken by the SPP code. 

The following notes are intended to clarify Figure 2-2. 

Alternate linkage data read in under the control of subroutine PR0B 
is read in in sequence (See Volume III, Section 2.11). 

Linkage data written by the various loss modules is positioned by sub- 
routine SETUP before being read by subroutine PICKUP. This insures 
that none of the linkage data is overwritten unless multiple cases are 
executed (in which case, the linkage data for the first case is destroyed). 
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3,0 Overlay ani Common Block Description 


Section 3,1 describes the overlay structure of the SPP while Section 3.2 
describes the common blocks used in the SPP code. 


3,1 Program Overlay Structure 

The SPP Program overlay structure is shown in Table 3-1 ai\d corresponds 
to the order of appearance of the subroutine write ups in Section 5, if read down- 
ward and then to the right. 

This overlay structure was developed under the restrictions of the CDC 
6000 series loader (SCOPE 3.3) which allows only three levels of overlay. Con- 
siderable savings in central memory storage can be achieved with a more exten- 
sive overlay structure is used. 
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BLOCK DATA 
ATMOS 



Table 3-1 Program Overlay Structure 
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3.2 Common Blocks 

The common blocks and the subroutines that are referred to by, are listed 
in Table 3-2. Tables 3-3 through 3-14 list the variables in some of the more im- 
portant labeled common blocks used in the SPP code. These tables are organized 
alphabetically by common block name. 

Appendix B of Reference 1 contains a dictionary of the variables appearing 
in labeled common blocks in the 0DE module. While there are somejninor vari- 
ances between Reference 1 and the SPP code, this dictionary should be very use- 
ful to anyone attempting to modify the ODE module. Unfortunately, similiar 
dictionaries are not available for the other modules. 


Table 3-2 Common Block References 


COMMON 

BLOCK 

REFERRED TO 

BY SUBROUTINE 

ADD 

BLKOAT 


CO'JTR* 


figi 


FIG* 


LOOPS 


PRISM? 


UNION 

amachr 

FAM 


TBLEDG 

ANMft 

PRINT 


BIKDAT 


DRIVER 


ODES 


UOKINP 


OUTPUT 


OUT 1 


PACK 


REACT 


RKTOUT 


SEARCH 

ARPRNT 

NZGEQP 


ODES 


ookinp 


PICKUP 
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COMMON 

BLOCK 


REFERRED TO 
BY SUBROUTINE 


BAUOH 

BALL 

CARD 

CDELMX 

cniF 

CIDAU 

CINT 

CK 

CQOATA 


SUMRY 

AVEBAL 

NZCF.Oh 

PICKUP 

SURRY 

BALLS 

SLKDAT 

LOOPS 

MAINJ 

PVS 

READIN 

READIN 

STORE 

DRIVER 

ODES 

SAVE 

START 

zetait 

BARCON 

BARPRO 

BARSET 

getpt 

QUITS 

READSI 

SEVAL 

START 

IAUX 
INT 
LESk 
M A I N 1 D 

ERROR 

READIN 

BARCON 

BARPRO 

barset 

EDI TCP 

GETPT 

QUITS 

READSI 

START 

zetait 





COMMON 

BLOCK 

REFERRED TO 

BY SUBROUTINE 

CODE 

(U.KDAT 


ERROR 


READIN 


STORE 

COOTOK 

ELU 


C0EF5X 

COF I IF 

COHCAS 


COMERS 

COMXP 

CQMy 


MA 1 N 1 D 

DDK 

ODKIhP 

CPHS 

SEARCH 

BARPRO 

PIU 

QUITS 

CONVRT 

OERIV 

DRIVER 

EF 

flu 

gtf 

IAUX 

INT 

maimd 

OPKINM 

OUTPUT 

PACK 

PRES 

PEAXIN 

SELECT 

STf 

T T APE 

DRIVER 

ODES 

CONvRr 

DERIV 

DRIVER 

EF 

*LU 

GTF 

IAUX 

InT 

MAIN1D 
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COMMON REFERRED TO 

BLOCK BY SUBROUTINE 


CONI 

COM2 

COM J 

COM4 

C0M5 

CPEfP 

CPRNT 

C3EVAL 

C3PRXC 


ODKINM 

OUTPUT 

PACK 

PRE3 

»EAXIN 

STP 

BLKDAT 

TAB 

TABIN 

BLKDAT 

MAINJ 

READJN 


BLKDAT 

MAINS 

READIN 

BLKDAT 

TAB 

TABIN 

BLKDAT 

DATAIN 

READIN 

TAB 

TABIN 

CPH3 

ROCKET 

DRIVER 

INT 

MAIN1D 

OOKiNP 

PACK 

BARPRO 

BARSET 

GETPT 

QUITS 

READSI 

SF.VAL 

START 

CONVRT 

OERIV 

If 

PACK 
RE AX IN 
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COMMON 

REFERRED TO 

BLOCK 

BY SUBROUTINE 


SELECT 

STOICC 

cutil cn^vRT 

CPH3 

DRIVER 

FLU 

IAUX 

HAIN1D 

MUK 

odes 

OOK 

UOKINP 

OUTPUT 

OUT 1 

PACK 

PICKUP 

PRES 

REACT 

REAXIN 

PKTUUT 

ROCKET 

SEARCH 

select 

SPLN 
STF 
SUHRY 
TT APE 

c M all o»ivm 

FL'J 

« a : m i n 

OHK INP 

PACK 

PRES 


OOLOUP 

DEPIV 


FF 


I AUX 


INT 


MAINJO 


flDKlNP t 


PRNTCK 

DOUBLE 

EDLHPM 


GAUSS 


MATRIX 


ODES 


OUT1 

E3GEX 

RARPPO 
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COMMON 

BLOCK 


REFERRED TO 
BY SUBROUTINE 


EDGCV 

m 

ERR 

EXTR^P 

FIG 

FnALBT 

MHTCOM 

1001 


REAOfI 

BARPfiQ 

GETPT 

REAOSJ 

PROBI.H 

SUHRV 

AD JK 

AGP 

CNTRt 

ERROR 

FCA'.C 

JArtES 

KPBPT 

N?HAIN 

PART II 

TD2P 

WALL 

«LPT 

DRIVER 

FIND 

PICKUP 

SURRY 


BLKDAT 

CONTRi 

CONTR2 

FIG 1 

FIG2 

PRISMl 

PRI3M2 

READIN 

STORE 

DERI V 
flu 

IAUX 

HAIN1D 

EQLBRM 

ODES 

AGP 

EISP 

ISP1D 









COMMON 

BLOCK 

REFERRED TO 

BY SUBROUTINE 

1N0X 

CPHS 


DRIVER 


EOLBRH 


FROZEN 


GAUSS 


hcalc 


MATRIX 


ODES 


OUTt 


REACT 


RKTOUT 


ROCKET 


SAVE 


SEARCH 

iNt'XX 

CPHS 


EQLBRM 

\ 

FROZEN 


GAUSS 


HCALC 


MATRIX 


OOFS 


REACT 


RKTOUT 


ROCKET 


SAVE 

W.RTS 

REA X IN 


SELECT 

I'JSt 

CPHS 


Driver 


EQLBR* 


FROZEN 


hcalc 


matrix 


ODES 


OUTi 


RKTOUT 

j 

i 

ROCKET 

1 

SAVE 


search 


SELECT 

K1NFU 

DRIVER 


frozen 


MAINJD 


ODES 


ODKINP 


OUTPUT 
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COMMON 

REFERRED TO 

BLOCK 

BY SUBROUTINE 

- 

OUT! 


PACK 


PRES 


REACT 


REAXIn 


RKTOUT 


ROCKET 


SELECT 


IKtM 

0ALIS 


MAlNj 


N2CEOM 


PICKUP- 


PROBLM 


8U M RV 

LKEQKN 

DRIVER 


ODES 


OOK 


PICKUP 


PROBLH 


RKTOUT 


ROCKET 


SURRY 

UKGEOH 

NZGEOM 


ODKINP 


SUHRV 


T02P 

LKMELT 

DERI V 


BAIN1D 

LKO^K 

PICKUP 


SUMRY 

LKT8L 

PROBIM 


READS! 

LKTD2P 

AGP 


NEGEOH 

• 

PICKUP 


PROBL* 


REAOSI 


T02P 

LOOP 

BALI 8 


BLKOAI 


C0NTR2 


ERROR 


PIG! 
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COMMON 

BLOCK 


REFERRED TO 
BY SUBROUTINE 








COMMON 

BLOCK 

REFERRED TO 

BY SUBROUTINE 


ROCKET 


Save 


SEARCH 

MISCX 

CPHS 


EQLBRM 


FROZEN 


HCALC 


MATRIX 


OOES 


REACT 


RKTOUT 


ROCKET 


SAVE 

MOUHTS 

OUT 1 


ROCKET 


SEARCH 

MUST 

OERIV 


FLU 

NAMD 

AGP 


EISP 


!8P1D 


TD2P 

NAME* 

A0CALC 


CCALC 


^CALC 


JAMES 


PAR TIL 


PCALC 


PROP 


TD2P 


TRACE 
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COMMON 

BLOCK 

REFERRED TO 

BY SUBROUTINE 

NAHfp 

ACOMP 


AGP 


AXI9PT 


CNTRL 


CON3T3 


EFN 


EISP 


I3P1D 


KPBPT 


ON tO 


PARTlL 


PRINT 


PTINT 


3UMPJ 


SUMP2 


TBUDG 


T02R 


TRACE 


wall 


«DG! 


WLPT 

NAMEr 

ACO*P 


ONED 


PARTIL 


TD2P 


TRACE 

NAMEI 

ADJK 


CNTRL 


KP9PT 


PRINT 


T02P 

NAHEL 

AGP 


ONED 


PARTI L 


PROP 


TD2P 


TRACE 


wDGI 
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COMMON 

BLOCK 

REFERRED TO 

BY SUBROUTINE 

NAHEh 

ACOHP 


AXI3PT 


CHECK 


CNTRL 


COASTS 


KPBPT 


ONID 


PARTIt 


PRINT 


PTINT 


8UMP1 


8UMP2 


T02P 


TRACE 


WALt 


WDGI 


WLPT 

NAmLN 

ACOMP 


AxISPT 


EEN 


kpbpt 


PTINT 


SUMPt 


SUMP2 


T02P 

NAMEP 

ACOMP 


ADJK 


AXI3PT 


CHECK 


CNTPL 


COASTS 


EEN 


KPBPT 


PRINT 


PTINT 


SUMP J 


8UHP2 


T02P 


Wl.PT 
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COMMON 

REFERRED TO 

BLOCK 

BY SUBROUTINE 

NAPL’Q 

ACOMP 

CNTRL 

COASTS 

ON tO 

PARTIL 

PRINT 

Toap 

TRACE 






NAHrft 

ONtO 

partil 

PROP 

TD2P 

TRACE 

NAMES 

ONED 

PARTIL 

PROP 

Toap 

TRACE 

WOGI 

NAMET 

AXI3PT 

KPBPT 

PTINT 

SUMP! 

SUMP? 

TD2P 

NIPT 

NAMEW 

QDKINP* 

Tl'ZP 


WALL 

NAME* 

ONED 

PARTIL 

PRINT 

PROP 

T8LEDG 

TD2P 

TRACE 

nahev 

ONED 

PARTIL 

PROP 

TD2P 

TRACF 

NAMEZ 

PARTIL 

Toap 

TRACE 






COMMON 

BLOCK 


REFERRED TO 
BY SUBROUTINE 


NAME! 

NAME2 

NAME3 

NAME'* 

NAMIN 

NIN 


ABC ALC 

CCALC 

rcAtc 

JAMES 

ONED 

PARTU 

PCALC 

TD2P 

ABCALC 

CCALC 

ecalc 

JAMES 
PC ALC 

ABCALC 

CCALC 

PCALC 

JAMES 

PARTIL 

PCALC 

PROP 

TRACE 

ABCALC 

CCALC 

MAIN1D ** 
ODKINP* 

BLKDAT 

ERROR 

REAOIN 

STORE 


ODECUT 


OUT 1 

RKTOUT 

ROCKET 
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COMMON 

BLOCK 

REFERRED TO 

BY SUBROUTINE 

PERP 

DRIVER 


EQLBRM 


PROZEN 


ODES 


OJT 1 


RXTOUT 


ROCKET 

PEPPX 

iqlbrm 


PROZEN 


ODES 


OUT 1 


RXTOUT 


ROCKET 

PMPAN 

BAR3ET 

POINTS 

DRtvER 


EQLBRM 


PROZEN 


HCALC 


MATRIX 


ODES 


OUTJ 


RKTOUT 


ROCKET 

POINT* 

EQLBRM 


PROZEN 


HCALC 


MATRIX 


ODES 


RKTOUT 


ROCKET 

PTABLE 

PLU 


I AUX 


MAINtD 


ODK 


OCKINP 


PRES 

PT3AVE 

DRIVER 


DOES 


RDLARO 


DRIVER 

PROBtM 
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COMMON 

BLOCK 

RDJRF.X 

RVACO* 

SAVE 

3AV00E 

SHGAMA 

SPFCE3 

SPECfX 


REFERRED TO 
BY SUBROUTINE 


DRIVER 
ODES 
PACK 
REAX JN 
SEARCH 

select 

STOICC 

INT 

MAINID 

ODK 

DDK I HP 
OUTPUT 

CONTRJ 

FIG2 

PRXSM1 

DRIVER 

ODES 

DERI V 
ODK 

ODK IMP 

CPHS 

EOLBRM 

P ROZEN 

HCAIC 

MATRIX 

ODES 

OUT 1 

RK TOUT 

ROCKET 

SAVE 

SEARCH 


CPHS 

EOLBRM 

PROZEN 

HCALC 

MATRIX 

ODES 

OUT l 

RK TOUT 

ROCKET 

SA^E 

SEARCH 






COMMON 

REFERRED TO 

BLOCK 

BY SUBROUTINE 

3PINP0 

CPH3 

DRIVER 

EQLBRR 

ODES 

ODKINP 

OUT 1 

RACK 

READS! 

RKTOUT 

ROCKET 

search 

SELECT 

8UMRY 

TD2P 

TBLDTi 

DRIVER 

IAUX 

MAIN JO 

OOKINP 

TBL0T2 

DRIVER 

IAUX 

HAINID 

TBLRC 

PICKUP 

SUMRV 

TBRSAV 

ODKINP 

Rt AX IN 

8ELECT 

T020UT 

DRIVER 

PRINT 

TBIEDG 

T02P 

TD2RE 

PICKUP 

SUHRV 

TfcQfcC 

ODES 

THI3P 

PICKUP 

SUHRY 

tmwuat 

NZGEOH 

PICKUP 

8‘JMRY 

TD2P 


3-21 





COMMON 

REFERRED TO 

BLOCK 

BY SUBROUTINE 

T0T3C 

DRIVER 

IAUX 

mainid 

TRICON 

DRIVER 

8UMRY 

TRAJEC 

TSTA01 

DRIVER 

IAUX 

MAINJD 

TWOPHZ 

CONVRT 

FLU 

MAINID 

ODK 

ODKINR 

OUTPUT 

pack 

VERSON 

sumry 

WHERE 

CONTR1 

CUNTR2 

FIG1 

EIG2 

PRTSM1 

PRISM* 

WROOK 

ODK 

OUTPUT 

wtelow 

BARPRO 

DIRECT 

READSI 

ZDEBUG 

DEM IV 

STF 

ZTR*N 

DRIVER 

IAUX 

MAINID 
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COMMON BLOCK: BALCOM 


FORTRAN 

NAME 

TYPE 

DIMENSION 

Description 

TIM 

L 

150 

Time, tp point from BAL module 

PTC 

R 

150 

total pressure at the end of the grain, 
P ci , (PSIA) 

RSTAR 

R 

150 

r* (ft) 

XLST 

R 

150 

L* (In) 

ETACS 

R 

150 

*0*1 

XMDOT 

R 

150 

rrij, mass flow rate (slugs/sec) 

EMIEMP 

R 

150 

not used 

PCI 

R 

150 

not used 

XMT0T 

R 

150 

not used 

FID 

R 

150 

not used 

PMAX 

R 

- 

n.aximum chambei pressure (PSIA) 

NBPT 

I 

- 

number of ballistics time points 

PAVE 

R 

- 

average P c (PSIA) 

RMIN 

R 

- 

smallest r* in the RSTAR array (ft) 

XLAVE 

R 

- 

average L* (in) 

ECAVE 

R 

- 

average r, c * 

XMAVE 

R 

- 

average m (slugs/sec) 



Table 3-3 Common Block: BALC0M 
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COMMON EEE 


F0RTRAN 

NAME 

TYPE 

DIMENSION 

DESCRIPH0N 


R 

15 

calculated ^ 


R 

15 

empirical rf 

ETAI 

R 

15 

Input tj 

ETAU 

R 

15 

value of v to be used 

IPETA 

I 

15 

pointer to value of ^to be used 


Table 3-4 Common Block EEF. 


i 
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C0MM0N BL0CK: 

: LKEQKN 

F0RTRAN 

NAME 


DIMENSION 

DESCRIPTION 

IPUISP 

I 

- 


If IPUISP =1; store € and ISP In 

AEQR and XIEQR 

IREQ 

I 

- 


If IREQ=1; restricted equilibrium op- 
tion; if 0, regular equilibrium 

AEQR 

R 

50 

< 

^ corresponding to I sp (XIEQR) for 
the restricted equilibrliim option 

XIEQR 

R 

50 



NEQR 

I 

- 


number of restricted equilibrium points 

SOLPT 

L 

- 


if TRUE solidification point proper- 
ties calculated 

FRZISP 

R 

- 


value of frozen ISP 



Table 3-5 Common Block LKEQKN 
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COMMON BLOCK: LK0DK 

FORTRAN 

NAME 

TYPE 

DIMENSION 

Description 

AKIN 

R 

40 

^ corresponding to XIK(I) 

XIK 

R 

40 

kinetic I calculated by ODK 

SD 

ETAKIN 

R 

- 

rtYihj at maximum area ratio specified 

irnOie $GEOM namelist 

NKPT 

I 

- 

number of points in the AKIN, XTK arrays 


Table 3-6 Common Block: LK0DK 


l 


I 

j 


« 


[ 
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C0MM0N BL0CK: LKTBL 


F0RTRAN 

TYPE 

DESCRIPTION 

NAME 



IF0DE 

I 

If IF0DE=1, ODE data Is available to the 
TBL module 

IFTD2P 

I 

If IFTD2P=1, TD2P data is available to 



the TD2P module 

INUN 

I 

logical unit on which the TD2P data Is 
available 


Table 3- 7 Common Block LKTBL 


1 



I 
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C0MM0N HL0CK: MAN 


J 

E 

l 


f 

I 


F0RTRAN 

DESCRIPH0N 

NAME 


0DF 

f if value = 0.0 not considered 

BAL 

1 if value =1.0 module by that name is to 

0DK 

J be executed 

TD2P 

^ if value = 2.0 module data was previously 

j generated 

TBL 

V if value = 3.0 loss to be read in 


Table 3-8 Common Block MAN 


| 


l 

I 
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C0MM0N HL0CK: MARKET 


F0RTRAN 

NAME 

TYPE 

DESCRIPTION 

LODE 

L 

If .TRUE, module 

ODE 

is to be executed 

LODK 

L 

if .TRUE, module 

ODK 

is to be executed 

LTD2P 

L 

if .TRUE, module 

TD2P 

is to be executed 

LTBL 

L 

if .TRUE, module 

TBL 

is to be executed 

LBAL 

L 

if .TRUE, module 

3AL 

is to be executed 


Table 3-9 Common Block MARKET 




I' 

l 
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C0MM0N BL0C: 

TBLRE 

FORTRAN 

NAME 

TYPE 

DIMENSION 

DESCRIPTION 

NTBL 

l 

- 

number of TBL points 

ABL 

R 

150 

A/A* from TBL 

DELISP 

R 

150 

A sp TBL 

DISP 

R 


^1 at the maximum area 

sp TBL 

ratio specified In the $GEOM 
namelist 


Table 3-10 Common Block: TBLRE 


i 
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COMMON BLOCK: TD2RE 

FORTRAN 

NAME 

TYPE 

DIMENSION 

Description 

A2D 

R 

150 

q corresponding to values of XI2D and E2D 

XI2D 

R 

150 

I ((.) from TD2P 

E2D 

R 

150 

V TD2P W 

NTDPT 

I 


number of points in the A2D, XI2D and E2D 
array 

ETD2P 

R 

- 

time average value of t)td2P 

XITD2P 

R 

- 

time averaged value of I 

sp TD2P 

ETD2PA 

R 

- 

7 j T D 2 P corr ected for throat erosion 

CD 

R 

- 

nozzle discharge coefficient as calculated 
by TD2P 

WD0T2D 

R 


nozzle mass flow rate as calculated by 

TD2P 


> 

i 

■ft _ ■ 
( 
t' 


i. 

f i 

! 


I 

I 


Table 3-11 Common Block: TD2RE 
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COMMON THISP 


F0RTRAN 

NAME 

TYPE 

DESCRIPTION 

ISP 

R 

Theoretical I at maximum area ratio 



specified In me $GEOM namelist 

CSTAR 

R 

Theoretical c* (ft/sec) as calculated by 



the 0DE module 


Table 3-12 Common Block THISP 
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COMMON HLOCK: THROAT 


FORTRAN 

NAME 

TYPE 

DESCRIPTION 

DRSDT 

R 

dr*/dt (ft/sec) 

TBIJRN 

R 

bum time (sec) 

NAR 

I 

number of points, at which nozzle erosion data is known 

FSUBM 

! R 

length of submergence/length of internal motor 

AEAS 

R 

nozzle entrance area/nozzle throat area 

AVELS 

R 

motor average L* (in) 

KN0S 

I 

nozzle type flag; =1 steel nozzle, =0 all others 

FCONP 

R 

, mole fraction of condensed phase (moles/100 grs of product) 













COMMON BLOCK: TRJCOM 

FORTRAN 

NAME 

TYPE 

DIMENSION 

Description 

ITJF 

I 

- 

trajectory flag, if 0 no trajectory input 

PAMBI 

R 

51 

ambient pressure (PSIA) 

TTRJ 

R 

51 

times corresponding to PAMBI 

NTJPT 

I 

- 

number of trajectory points input 


Table 3-14 Common Block TRJCOM 


r 

f. 




I 
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4.0 Program Fi les 


The SPP program uses nine files to read, write, punch and temporarily save 
data. Six of these files are always referred to by variable names and hence, may 
be changed to suit the restrictions of the computer system being used. These 
variable names are set in Subroutine DRIVER for units 25 through 29 and by input 
data for unit 3 in subroutine NZGE0M. 

Logical unit 8 is equated to the system punch file for CDC 6000 series 
computers and it is rewound. Therefore, care must be taken in assigning this 
file on machines which do not allow rewinds on the system punch file (e.g., Univac 
1108). This linkage data written out on this unit is described in Volume III, Section 
2.11 and will not be repeated here. 

Table 4-1 lists the files used by the SPP code and the FORTRAN variable 
names used to reference these files. The device types listed correspond with the 
current usage on CDC 6000 series computers. The mnemoics used in Table 4-1 
are listed below. 

T - Tape files 

MS - Mass Storage (perm) files 

TMS - Temporary Mass Storage (scratch) files 
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FORTRAN 

VARIABLE 

NAME 

Loc^cal 

TJnit 

Device Type 

Description 

INP 

3 

T, MS, TMS 

alternate Input file for SPP linkage 
data 

- 

5 

card reader 

data input sLream 

- 

6 

line printer 

printed output stream 

- 

8 

punch 

module linkage data output file 

JANAF 

25 

T, MS, TMS 

thermodynamic data file 

KREAX 

2f> 

TMS 

temporary file used to store reaction 
rate data 

KSTF 

27 

TMS 

temporary file used to store a short 
version of the thermodynamic data 

IRREAD 

28 

TMS 

temporary file used for multiple cases 
in tne 0DK module. 

ITSTAB 

29 

TMS 

the data written on this file is not 
currently used 


Table 4-1 Files U3ed By The SPP Code 


i 
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5.0 PROGRAM SUBROUTINE DESCRIPTIONS 


This section of Volume II contains a subroutine by subroutine description 
of the SPP computer code. Dictionaries of some of the more important variables 
which appear in labeled common areas of the program appear in Section 3. 

Section 5 is organized in accordance with the overlay structure of the SPP 
code and the user is advised to refer to Section 3 of this document for the location 
of individual subroutines and common blocks and their cross references. 


5-1 


5.1 Link 00 Subroutines 


5.1,1 Program Driver 

This is the main routine for the SPP code and as such provides the overlay 
control , defines most of the upper level labeled common blocks, and initializes 
most of the control variables used in the code. On CDC 6000 computers this 
routine also defines all of the files and buffer sizes which are used in this code. 

The DRIVER routine is essentially driven by directive cards which are read 
in on logical unit 5. The directive cards which are recognized by this subprogram 
are listed below in Table 5.1-1 along with the action taken by the DRIVER routine. 


Table 5.1-1 Data Input Directives 


Directive card 

Action Taken 

THERM 0 

Call in overlay 1,0 which processes ther- 
modynamic data. 

TITLE 

Stores the case title for subsequent print 
out. 

PR0BLEM 

Cells in the PR0BLM subroutine which de- 
termines which loss modules are to be ex- 
ecuted. This directive and subsequent 
data initiate the execution of the various 
loss modules. 

LOW T CPHS 

Calls in subroutine LTCPHS which processes 
low temperature thermodynamic data. 

GE0METRY 

Calls in subroutine GE0M yvhich processes 
the nozzle geometry data. 

TRAJECT0RY 

Call in subroutine TRAJ which processes 
the trajectory and/or ambient pressure 
history of the motor. 


In addition to processing the directive cards and calling in the selected 
loss modules, the main program also hi idles the communication of the linkage 
data between modules and calls the summary data output routine SUMRY. 
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5. 1 . 2 Subroutine BL0CK DATA 


BL0CK DATA contains atomic data stored in AT0M(i,j) and many of the vari- 
ables used with the variable format, FMT. The AT0M variables are defined in 
appendix B, Reference 1. The format variables are stored in the common labeled 
0UPT and are described here. 

A variable format was used so that one format, FMT, could be used in the 
final output with changes in the number of decimal places according to the sizes 
of the numbers. The format is used to print a label and from 1 to 13 associated 
numbers. The labels contain 14 alphameric characters stored in four words and 
printed with 3A4, A2. The numbers are all printed in a field of 9. FMT is initially 
set in HL0CK DATA as follows: 


FMT 

(i) 

(2) 

(3) 

(4) 

(5) 

(6) 

(7) 

(8) (9) 

do) 

(11) 


(1H 

,3A4 

,A2, 

F9. 

0 , 

F9. 

0 , 

F9. 0.. 

F9. 

0 , 

FMT 

(12) 

(13) 

(14) 

(15) 

(16) 

(17) 

(18) 

(19) 

(20) 

(21) 


F9. 

0 , 

F9. 

0 , 

F9. 

0 , 

F9. 

0 , 

F9. 

0 , 

FMT 

(22) 

(23) 

(24) 

(25) 

(26) 

(27) 

(28) 

(29) 

(30) 



F9. 

0 , 

F9. 

0 , 

F9. 

0 , 

F9. 

0 

) . 



where the spaces are stored as blanks. 


Some variables set in BL0CK DATA to modify FMT are as follows: 
Variable: FO FI F2 F3 F4 F5 FB FMT13 FMT9X FMT19 

Storage: 0, 1, 3, 4, 5, 13, 9X, 19, 

The following is a list of variables used as labels and printed with 3A4, 

A2 in FMT: 


Variable 

Stored label 

FP 

P, ATM 

FT 

T, DEG K 

FH 

H, CAL G 

FS 

S, CAL/(G) (K) 

FM 

M, MOL WT 

FV 

(DLV/DLP) T 

FD 

(DLV/DLT) P 

FC 

CP, CAL/ (G) (K) 

FG 

GAMMA (S) 

FL 

SONVEL, M/SEC 

FR1 

PC/P 

FC1 

CF 

FN 

MACH NUMBER 

FR 

CSTAR, FT/SEC 

FI 

ISP, LB-SEC/LB 

FA 

IVAC, lb-secab 

FA1 , FA2 

AE/AT 
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5.1.3 Suoroutine ATM0S 


Soubroutine ATM0S is used by subroutine TRAJEC (see Section 5.1.16) to 
calculate the ambient pressure history during a motor firing if the ambient condi- 
tion time history is given as a function of altitude. 

This subprogram calculates, given the geometric altitude, the ambient temp- 
erature PR), pressure (psia), sound speed (ft/sec), density (slugs/ft 3 ), and pres- 
sure gradient (lb/ft 3 ). Subroutine ATM03 assumes that the atmosphere consists 
of a gas of constant molecular weight on a spherical earth whose atmosphere is 
in a hydrostatic balance. The temperature profiles used in the hydrostatic balance 
are those given in the 1962 ARDC standard atmosphere tables. 

These assumptions yield data which agree with that standard atmosphere 
table to within 1% up to altitudes of 700 kilometers (12,296,585 ft.). 


5.1.4 Subroutine AVEBAL 

Subroutine AVE8AL is called by subroutine PICKUP and calculates from the 

linkage data supplied by the internal ballistics calculation the time averaged values 

of chamber pressure, L*, throat radius, rj *, and mass flow rate. Also printed out 

c 

are the maximum and minimum values of the above quantities along with their time 
histories. The throat erosion rate in inches/second is also calculated and printed 
out. 
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5.1.5 Subroutine ETACFX 


This subroutine calculates the nozzle performance losses using empirical 
correlations. The losses calculated in this routine are in terms of per-cent nozzle 
thrust coefficient lo ss. Subroutine SUMRY converts these losses into fractional 
I g p efficiencies using the equation 

100 ^ 

*ISV 1 = Too *»c* 5,1_1 

where ^is the nozzle Cp loss as calculated below 

Care should be taken when using the above equation that c* efficiency, 

n *, is not used more than once to convert C r efficiencies to I efficiencies. 

"C 1 * r Sp 

The losses considered and equations used by this routine are given below. 

A more complete description of these losses is given in Voi. I, Section 4. 


Nozzle divergence loss is calculated as: 



5.1-2 


where <y=half angle of conical nozzle, included angle of contoured nozzle 
0Ex =ex it angle of contoured nozzle 

The loss due to finite rate chemical kinetics is calculated as: 


^KIN = 33,3 


theoretical frozen Isp ,200 \ 
theoretical shifting Isp * P 


where P=nczzle chamber pressure in psia or 200, whichever is greater. 


The boundary layer loss is expressed in the form of a time-dependent heat 
transfer expression, with a correction term for nozzle expansion ratio: 


, 0.8 


^BL " C 1 


D 


0.2 


1 + 2 exp (-C 2 P 0,8 t/D°' 2 ) 


1 + 0.016U- 9) 


5.1-4 


where P = pressure, psi 

= throat diameter, in 
t - time, sec 
e - expansion ratio 
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The coefficients and C 2 are determined as Is shown below. The nozzle 
type is input by the user as KN0Z in the $GE0M namelist. 

Ordinary Nozzles: = 0.00365 

C 2 = 0.000937 

Steel Nozzles: = 0.00506 


The nozzle two phase flow loss term is calculated as: 

? C 4Dp C5 

^TP U 3 p 0.15 f 0.08 D C g 

where £ = mole fraction of condensed phase, moles/100 gm 

Dp= particle size, u 
P = pressure, psi 
<r •= expansion ratio 
= throat diameter, in. 
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The coefficients are determined from the following: 
If £ ^ 0.09, C 4 = 0.5 



with D p =• 0.454 P 1//3 £ 1//3 [l - exp (-0.004L*) J [ 1 + 
and L* = motor characteristic length, in. 

The motor submergence loss ‘is calculated by: 

pi, 0.8 qO • 4 

"SUB = 0 * 0684 ( A^ XT 

U t 

where 

P = pressure, psi 

£ = mole fraction of condensed phase, moles/100 gm 

A* = nozzle entrance area/nozzle throat area 
S = length of submergence/length of internal motor 
D t = throat diameter, in. 


i 

f 

t 


0.045D t J 5.1-6 


5.1-7 
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5,1,6 Subroutine FIND 


This subroutine locates the indes, I, in a table such that X(I) <X <X(I+1). 
On option, when the variable IEXTP in labeled common IXTRAP^O, subroutine FIND 
allows for extrapolation in tables. See Appendix A for conversion aides. 


5,1.7 Subroutine LINK1 A 

Subroutine LINKlA's only function is to call subroutine PR0BLM. It is 
provided to allow addition overlay capability on computers with explicid overlay 
calling features. 


5.1.8 Subroutine L7CPHS 

This subroutine processes the low temperature C , H, S Thermodynamic 
Data extension input as described in detail in Volume III, Section 2.2.1. The 
low temperature values of C p , H, and S are at present only used in the 0DK 
module. 
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5,1,9 Subroutine NZGE0M 

Subroutine NZGE0M Is called by program DRIVER when a GE0METRY directive 
is encountered In the card input stream. This routine uses the $GE0M Namelist 
input to define the nozzle geometry. Data which are read under this Namelist 
are stored in labeled common blocks LKGE0M, LKBM, ARPRNT, BALC0M, LKTD2P, 
and THR0AT. 

The variables in the above common blocks are described in Section 3,2, 
However, for completeness, the variable names whose values are determined in 
subroutine NZGE0M are listed below in Table 5,1-1, 


C0MM0N HL0CK 

VARIABLE NAME 

ARPRNT 

ASUB 


AS UP 


NASUB 


NAS UP 

BALC0M ~~ 

RSTAR 


TIM 

LKBM 

NN0Z " 


R 

LKGE0M 

GMVAR 


IGM 


IWF 


NWS 


RS 


ZS 

LKTD2P 

RSTARM 


XLS 

THR0AT 

DRSDT 


FSUBM 


NAR 


TBURN 


Table 5,1-1 Common Variables Set in NZGE0M 

The variables shown in the above table are communicated downward only 
in the SPP code. That is, values set in this routine may be changed in indivi- 
dual calculational mcdules, however, the values which are changed are not 
transmitted to other calculation modules. Instead the initial values set via the 
$GE0M Namelist are used. 
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5.1,10 Subroutine PICKUP 


Subroutine PICKUP processes all of the linkage data between the five basic 
computational modules which are transmitted on logical unit 8 except for the 
boundary layer edge conditions generated by the TD2P module for the TBL module. 

This routine also performs unit conversion on the data read in, checks for 
overflow of tables, and prints out a summary of the results for each of the modules 
which have been executed. 


5.1.11 Subroutine FR0BLM 

Subroutine PR0BLM is executed by program DRIVER when the PROBLEM 
directive is encountered in the card input stream. 

This routine sets flags and determines via user input in the $PR0B Namelist 
which of the five basic computational loss modules are to be executed. If the 
option to read previously generated data is selected, subroutine PICKUP is called 
by PR0BLM to read these data. 


5.1.12 Subroutine SETUP 

Subroutine SETUP is used to position the linkage data file, logical unit 8, 
for processing by subroutine PICKUP. This routine reads through the linkage data 
file until the appropriate calculation module name is found, then it backspaces 
one logical record so that the loss module name is the first card image read by 
subroutine PICKUP. 
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5.1.13 Subroutine SPLN 


Performs either cubic or linear Interpolation between two given points. 


Cubic interpolation for a function and Its first two derivatives Is performed 
as described below: Given function values y n and y fi+1 and first derivative values 
y' n and y^ +1 at x R and x n+1 , this subroutine evaluates y(x), y'(x), and y"(x) for 
x n <:x<x n+1 using: 

y = A(x - x J 3 + B(x - x ) s + C(x - x ) + D 
n n n 



V- (y; H - v') A 

where - 

A = F • [^n+l + ^ h -*J 

B = ”F * [ ' y n+l + 2y n )h ' 3k ] 

c y ; 

D= y n 
h = Vl • x n 


Linear interpolation for a function and its first two derivatives is performed 
as described below: 



V' 


Vl I y n 

x n+l " X n 


y" = 0.0 
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5.1,14 Subroutine SUMRY 


This subroutine is called by program DRIVER (Section 5.1.1) and prints out 
a summary report on the results of each of the five basic computational modules 
(if executed) at the maximum area ratio requested in the $GE0M namelist. If 
the GE0METRY directive was not input to the SPP code the SUMRY routine assumes 
that one or more of rhe calculation modules were executed independently and 
that no summary of results is desired. 

There are three basic sources for each of the losses considered by the SPP 
code. These sources are due to one of the following items shown in the table below. 


1. User Input 

2. Execution of a loss module 

3. Calculation of an empirical correlation 


Table 5.1-3 Sources of Losses Considered 

User input of loss efficiencies along with selection of which of the loss 
modules are to be executed is determined by input to subroutine PR0BLM (Section 
5.1.11). Empirical losses are calculated by subroutine ETACFX (Section 5.1.5) 
which is called by SUMRY. While all of the above losses are printed out in this 
routine, the selection of which source of the loss is to be used in determining the 
total delivered perfoimance of the solid propellant rocket motor is given in the 
order appearing in Table 5.1-3. That is, if the user of the SPP code inputs a loss, 
that is the loss efficiency that is used in calculating the final delivered perfor- 
mance. Hence, while user input overrides values calculated by the analytical 
loss modules, the values calculated by the loss modules takes precedence over the 
values determined by the empirical correlations. 

Using the selection criteria described above, subroutine SUMRY then prints 

out the product of efficiencies along with the decrement of performance due to a 

turbulent boundary layer. Hence, if no input efficiencies were input, the value 

of the delivered performance due to input losses would then be just the theoretical 

I calculated by the 0DE module, 
sp 




In order to determine the motor delivered performance to ambient conditions, 
trapezodial rule integration is used to calculate the average ambient pressure and 
average expansion ratio if this information is available. That is, the average 
performance delivered to ambient conditions is calculated by 

_ n 

I =1 , fj. i, - - aI 5.1-8 

sp D sp th 1=1 ‘ l ' k sp TBL sp A 

where I = I calculated by 0DE 

sp th sp VAC 

Yf . = fractional loss efficiency determined via the selection rule, 

k, for the I th loss considered 

A I = loss decrement due to the formation of a turbulent boundary 
sp TLL layer 

A I = loss decrement due to the motor delivering thrust to an am- 

sp A blent pressure. This value is calculated as (P . /?) 

+ // r* \ amD c 

x j xcv^xey 

If both the ballistics and the two dimensional -two phase flow modules link- 
age data are available, then subroutine SUMRY calculates the delivered perfor- 
mance throughout the motor firing for each time point specified in the trajectory 
and ballistics input data. The following formulas are used to compute the values 


printed out by this routine. 

n CE (t) = max (1.0, C D * tr t .*(t)) 

sp TD2P^ 

Isp VAC (t) = Isp th * f, ® 2pW ' 4 sp TD2p (t=0) ) * ’ 7CE<t) ‘ ^ ' ^TBL 


F VAC (t) 

= I * m(t) 

sp VAC 

F amb (t) 

= P amb^ - A exlt 

Isp amb (y 

= F a , b (t)/m(t) 

f del 

= F VAC (t) " F amb (t) 

Isp DEL 

Isp VAC Isp amb 

I 

= J F DEL ^ dt 

0 

m 

= J* mdt 


0 
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where 


m 


n c* 


amb 


^SUB 

%IN 




sp bl 


^TD2P 

Sp TD2P 

A .* 
exit 

*amb 

I 


m 

I 

sp DEL 

F 

DEL 

I 

sp amb 


mass flow rate (from Ballistics calculation) 
c* efficiency (from Ballistics calculation) 
ambient pressure (from input tables) 

submergence loss efficiency (selected as described above) 
kinetic loss efficiency (selected as described above) 
boundary layer loss decrement (selected as described above) 

2D-2phase flow loss efficiency (from TD2P as a function of *(t)) 

2D-2phase I (from TD2P as a function of ^ (t)) 

sp 

no?zle exit area (from input) 

thrust decrement due to ambient pressure 

total impulse to this time 
total mass expended to this time 
delivered specific impulse 

delivered thrust 

I decrement due to ambient pressure 
sp 
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5.1,15 Subroutine TIMERX 


This subroutine Is provided to localize printing of computer execution times 
during the calculations. This subroutine should be modified at each local instal- 
lation to provide proper interface with local system timing routines. This subroutine 
calls a CDC system subroutine named SEC0ND to obtain the current run execution 
time in seconds. 


5.1.16 Subroutine TRAJEC 

Subroutine TRA/EC is used to calculate ambient pressure during a motor firing. 
A table of ambient pressure or altitude is input to this routine via the Namelist 
name $TRAJ as a function of time. Delivered performance to the specified amiJent 
conditions is then calculated in subroutine SUMRY and printed. If an altitude 
versus time table is input to subroutine TRAJEC, then subroutine ATM0S is used to 
calculate the ambient pressure at which the performance of the motor is delivered. 
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5.2 Link 10 Subroutines 


5.2.1 Subroutine LINK10 

This subroutine consists only of a call to subroutine TTAPE and is used 
only to facilitate conversion of the program overlay structure to other computers. 


5.2.2 Subroutine TTAPE 

When a THERM# directive card is read by the main program this subroutine 
is called to generate a master Thermodynamic Data tape (Logical Unit 25). The 
input Thermodynamic Data is in curve fit form and is identical to that required for 
the #DE computer program described in NASA SP-273, Reference 1. The format 
for this curve fit date is described in the Volume III, Section 2.2. This routine 
does not update the master Thermodynamic Data tape, but instead creates a new 
one. Hence, all of the thermodynamic must be read in, old and new, if the user 
wishes to add new species to this file. 
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5.3 Link 20 Subroutines 

Overlay Link 20 contains the 0DE module which calculates the theoretical 

or reference I used in the methodology incorporated in the SPP code, 
sp 

The main body of subprograms which are contained in this overlay are from 
the NASA-Lewis Equilibrium Program (Ref.l) as modified for rocket performance 
calculations (Ref. 3). The brief analytical descriptions of the equilibrium calcu- 
lations contained herein were condensed from Reference 1. For a complete des- 
cription of these calculations see Reference 1 (NASA SP-273). 


5.3.1 Subroutine LINK20 

This subroutine consists only of a call to subroutine 0DES and is used only 
to facilitate conversion of the program overlay structure to other computers. 
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man 








5.3.2 Subroutine CPHS 


Cp° 


This subroutine evaluates the thermodynamic functions 


H X 


S °T 


RT ' R 

from curve fit coefficients. Two sets of coefficients are used for two adjacent 
temperature Intervals. The functions evaluated are presented below: 

Cp\ 


~ 1T 

h° t 

~kT 

S T 
R 

G °T 

TT 


a, + a 0 T + a.T 2 + a ,T 3 + a.T 4 
1 2 3 4 5 


a i + 


a 2 T 


a 3 T " 


a 1 InT + a 2 T + 


a 3 T " 


a.r 


a.T 3 

4 


a 5 T4 


T 


a.T 4 

+ " i 4 _ + a 7 


H °T 

RT 


S° 


R 


5.3.3 Subroutine DETON 

This is a dummy subroutine used to replace subroutine DET0N of Reference 

1. 


5.3 .4 Subroutine EFMT 

Subroutine EFMT (E-format) writes statements in a special exponent form. 
This form is similar to the standard FORTRAN E-format, but the letter E and some 
of the spaces have been removed for compactness. This routine is not used in 
the SPP code at the present time. 
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5,3,5 Subroutine EQLBRM 


EQLBRM is the control routine for the equilibrium module which calculates 
equilibrium compositions and thermodynamic properties for a particular point. A 
free-energy minimization technique Is used. This routine permits calculations 
such as (1) chemical equilibrium for assigned thermodynamic states (T, P) , (H,P), 
(S,P), (T,V), (U,V), or (S,V), (2) theoretical rocket performance for both equili- 
birum and frozen compositions during expansion, (3) incident and reflected shock 
properties, and (4) Chapman- Jouguet detonation properties. The program considers 
condensed species as well as gaseous species. A detailed description of the 
equations and computer program for computations involving chemical equilibria 
in complex systems is given in Reference 1. Figures 4(a) through 4(c) of Reference 
1 give a complete flow diagram for this subroutine. 


5,3,6 Subroutine FR0ZEN 

Subroutine FR0ZEN is called from subroutine R0CKET to calculate the 
temperature and thermodynamic properties for the following assigned conditions: 

1. Composition frozen at combustion conditions. 

2. At assigned exit pressure. 

3. At assigned entropy equal to the entropy at combustion conditions. 

The iteration procedure used for obtaining the exit temperature is discussed in the 
section Procedure for Obtaining Frozen Rocket Performance (p. 40, Reference 1). 

This subroutine has been modified to calculate the strictly frozen performance 
at well below the solidification point of any condensed phases which are present. 
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5.3.7 Subroutine GAUSS 


Subroutine GAUSS Is used to solve the set of simultaneous linear Iteration 
equations constructed by subroutine MATRIX. The solution Is effected by perform- 
ing a Gauss reduction using a modified pivot technique. In this modified pivot 
technique only rows are interchanged. The row to be used for the elimination of 
a variable is selected on the basis that the largest of its elements, after divi- 
sion by the leading element, must be smaller than the largest element of the other 
rows after division by their leading elements. 

The solution vector is stored In X(k), In the event of a singularity, IMAT 
(which is equal to the number of rows) is set equal to IMAT - 1. IMAT is tested 
later in subroutine EOLBRM. 


5 . 3 . 8 Subroutine HCALC 


The purpose of HCALC is to calculate thermodynamic properties for reactants 
under certain circumstances. HCALC is called from entry NEW0F of SAVE. 


HCALC is called from NEW0F when CALCH is set true . CALCH is set true 
in the main program when zeros have been punched in card columns 37 and 38 on 
one or more REACTANTS cards. The zeros are a code indicating that the enthalpy 
(or internal energy for UV problems) for the reactant should be calculated from the 
THERM0 data at the temperature punched on the card. This temperature has been 
stored in the RTEMP array. CPHS is called to calculate the enthalpy. The value 
is stored in the ENTK array and printed in the final tables. 

The properties calculated in subroutine HCALC, their FORTRAN symbols, 
and the conditior for which they are used are as follows: 


Property 

F0RTRAN Symbol 

When used 

H(K)T 

HPP(K) 

if cc 37 and 38 are 

V R 

HSUBO 

input as 00 on reactants card 


5.3.9 Subroutine MATCH 

This subroutine is called by R0CKET to supply subroutine MUK with a vector 
of internal sequence m.mbers which point to the appropriate Lennard-Jones parameters 
used by MUK to calculate transport properties. The input to MATCH is the species 
name (from the SUB array) and the output corresponds to the index numbers in the 
table in the MUK write up. 
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5.3.10 Subroutine MATRIX 


This subroutine- sets up the matrices corresponding to tables I through IV 
of Reference 1. The assigned thermodynamic state being set up (tables I and II) 
is specified by the following codes: 


Assigned 

thermodynamic 

state 

Codes 

TP 

TP = .TRUE. VOL = .FALSE. CONVG = .FALSE. 

HP 

HP= .TRUE. VOL = .FALSE. CONVG = .FALSE. 

SP 

SP = .TRUE. VOL = .FALSE. CONVG = .FALSE. 

TV 

TP = .TRUE. VOL = .TRUE. CONVG = .FALSE. 

UV 

HP= .TRUE. VOL = .TRUE. CONVG = .FALSE. 

SV 

SP = .TRUE. VOL = .TRUE. CONVG = .FALSE. 


After convergence of any of the previous six problems, setup of the deriva- 
tive matrices (tables III and IV) is specified by the following codes: 


Derivative 

Codes 

DLVTP 

DLVPT 

1 1 


Note: Only the Rocket option is available in the SPP code. 
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5.3.11 Subroutine MUK (CP, XW, T, C, IT, N, XM, XMT, TK, PR) 

This routine calculates the viscosity, thermal conductivity and Prandtl number 
for a gas composed of a mixture of species. 

CALLING SEQUENC E; 


CP 

is an array of N specific heats of the 
species (ftVsec 2 6 R) 

(INPUT) 

XW 

is an array of N molecular weights of 
the species (slug/slug mole) 

(INPUT) 

T 

is the temperature of the gas fR) 

(INPUT) 

C 

is an array of N mass fractions of the 
species 

(INPUT) 

IT 

is an array of N indices of the species 
into a table of collision diameter (cr) 
and energy of attraction (e/k) 

(INPUT) 

N 

is the number of species in the gas 

(INPUT) 

XM 

is an array of N viscosities of the 
species, (Ibf • sec/ft 2 ) 

(0UTPUT) 

XMT 

is the total viscosity of the gas, 

(lbf • sec/ft 2 ) 

(0UTPUT) 

TK 

is the thermal conductivity of the gas, 
(ft.lbf/ft 2 sec PR/ft)) 

(0UTPUT) 

PR 

is the Prandtl number of the gas, 

(0UTPUT) 


Method: 

N 

C ^ E C. C Specific Heat 

p i=i 1 p i 


M 


w 


-fir 

2 

i=i 




Molecular Weight 
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M 


» C. . 


w 


i * M 


w, 


*1 


T/1.8 

urn; 


fl = table (T*) 


table of n vs T* 


= table (i) 


«\/ M T 

4.15822x10" V w i 

/ M \ — 1/2 
1 ll+ Wl 


n r / n x. \ 
= Si Ui ( 1+ 5i x ij 

L X tfi / 


table of a and c/k 
vs, individual 
species 


». \ ^ i \ ' 1/1 



J 


-1 


W 


W, 


viscosity of the gas 


K, 


K 


M, R 


M 


(.45 + 1.32 


Pi 


w x 


rnr) 


) 


/ N 

!■ K £ f 1 + lo065 S 


N 

r. 

1=1 


j.) 

i i x i 1 


-i 


j=i i 

i^i 


thermal conductivity 


K 


Prandtl number 
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Equations for and ]U i are from Reference 7. The values of the collision 
integral are from table 2 of Appendix B in Reference 8. The relations used 
to calculate /land K of the mixture are from Reference 9 . 

'Also from Reference 7 are the values of the collision diameters, a, and energy 
of attraction, e/k. 

Table 1 correlates the chemical name of the species to the internal number 
assigned to it by the subroutine. Also included in Table 1 is the key-punch 
name assigned to 'die species, since lower-case letters and subscripts are 
non-standard features in most computer configurations. 
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Table 1 . SPECIES NAMES AND IDENTIFIERS 


Species 

Number 

Chemical 

Name 

Key Punch 
Name 

1 

A1 

AL 

2 

A1C1 

ALCL 

3 

AlCla 

ALCL3 

4 

A1F 

ALF 

5 

A1F 3 

ALF3 

6 

AIN 

ALN 

•/ 

A10 

ALO 

3 

A1S 

ALS 

9 

Ala 

AL2 

10 

Air 

AIR 

11 

Ar 

AR 

12 

AsH 3 

ASH3 

13 

B 

B 

14 

BBr 3 

BBR3 

15 

BC1 

BCL 

10 

BC1 2 

BCL2 

17 

BC1 3 

BCL 3 

16 

BF 

BF 

19 

bf 2 

BF2 

20 

bf 3 

BF3 

21 

bi 3 

BI3 

22 

BO 

BO 

23 

b(och 3 ) 3 

B(OCH3)3 



Species 

Number 

Chemical 

Name 

Key Punch 

Name 

35 

Br 

BR 

36 

BrF 

BRF 

37 

BrF 3 

BRF3 

38 

BrO 

BRO 

39 

Br 2 

BR2 

40 

C 

C 

41 

CBrF 3 

CBRF3 

42 

CBr4 

CBR4 

43 

CC1 

CCL 

44 

CC1F 3 

CCLF3 

45 

CC1 2 

CCL2 

46 

CC 12 F 2 

CCL2F2 

47 

CC1 3 

CCL3 

48 

CC 13 F 

CCL3F 

49 

ecu 

CCL4 

40 

CF 

CF 

51 

cf 2 

CF2 

52 

CF:j 

CF3 

53 

cf 4 

CF4 

54 

CH 

CH 

55 

CHBrCIF 

CHBRCIF 

56 

CHBrCU 

CHBRCL2 

57 

CHBr3 

CHBR3 

58 

CHC1F 2 

CHCLF2 

59 

CHClg 

CHCL3 

GO 

chf 3 

CHF3 

Cl 

CHpBrCl 

CH2BRCL 

62 

CHnCIF 

CH2CLF 

63 

CH S CL«, 

CH2CL2 

64 

CH 2 F 2 

CH2F2 

G5 

ch 3 i 2 

CH2I2 

66 

CH 3 Br 

CH3BR 

67 

CH 3 CI 

CH3CL 
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Chemical 

Name 


Species 

Number 

58 

69 

70 

71 

72 

73 

74 

75 

76 

77 

78 

79 

80 
81 
82 
03 

84 

85 

86 
37 
88 

89 

90 

91 

92 

93 

94 

95 
36 

97 

98 

99 


CH 3 F 

CH3I 

CH3OH 

CH4 

CN 

CO 

cos 

co 2 

CP 

cs 

CSa 

Cs 

QjHa 

C S H 

CsHq 

CjHgCl 

CsHgOH 

C 3 N 2 

CH3OCH3 

ch 3 chch 3 

CH3CCH 

cyclo-CsHg 

C3He 

n-CsH^OH 

CH3COCH3 

CH3COOCH3 

n-C 4 H 10 

iso-C 4 H 10 

CsHsOCbHs 

CHsCOOCsHg 

n-C 5 H 13 

C(CH 3 ) 4 


Key Punch 
Name 

CH3F 

CH3I 

CH30H 

CH4 

CN 

CO 

COS 

C02 

CP 

CS 

CS2 

C2 

C2H2 

C2H4 

C2H6 

C2H5CL 

C2H50H 

C2N2 

CH30CH3 

CH2CHCH3 

CH3CCH 

CYCLO-C3H6 

C3H8 

N-C3H70H 

CH3COCH3 

CK3C00CH3 

N-C4H10 

ISO-C4H10 

C2H50C2H5 

CH3COOC2H5 

N-C5H12 

C(CH3)4 



Species 

Chemical 

Key Punch 

N amber 

Name 

Name 

101 

C 6 H 12 

C6H12 

102 

n-CeHx* 

N-C6H14 

103 

Cd 

CD 

104 

Cl 

CL 

105 

C1CN 

CLCN 

106 

C1F 

CLF 

107 

C1F 3 

CLF3 

108 

CIO 

CLO 

109 

Cl 2 

CL 2 

110 

F 

F 

111 

FCN 

FCN 

112 

f 2 

F2 

113 

H 

H 

114 

HBr 

HBR 

115 

HCN 

HCN 

116 

HC1 

HCL 

117 

HF 

HF 

118 

HI 

HI 

119 

HS 

HS 

120 

h 2 

K2 

121 

h 2 o 

H2G 

122 

h 2 o 2 

H202 

123 

1 .-sb 

H23 

124 

He 

HE 

125 

Hg 

HG 

126 

Hg Br 2 

HGBR2 

127 

HgCle 

HGCL2 

128 

Ugh 

HGI2 

129 

I 

I 

130 

ICI 

ICL 

131 

h 

12 

132 

Kr 

KR 

133 

Li 

LI 

134 

LiBr 

LIBR 
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Species 

Chemical 

Key Punch 

Number 

Name 

Name 

135 

LiCN 

LICN 

136 

LiH 

LI CL 

137 

LiF 

LIF 

138 

Lil 

LII 

139 

LiO 

LIO 

140 

\j1l2 

LI2 

141 

LijjO 

1/20 

142 

Mg 

MG 

143 

MgCl 

MGCL 

144 

MgCla 

MGCL2 

145 

MgF 

MGF 

146 

MgF 2 

MGF2 

147 

Mg 2 

MG2 

148 

N 

N 

149 

NFa 

NF3 

150 

NH 

NH 

151 

nh 3 

NH3 

152 

NO 

NO 

153 

NOC1 

NOCL 

154 

N 

2 

N2 

155 

n 2 o 

N20 

156 

Na 

NA 

157 

NaBr 

NABR 

158 

NaCN 

NACN 

159 

NciJl 

NACL 

160 

NaF 

NAF 

161 

Nal 

NAI 

162 

NaO 

NAO 

163 

NaOH 

NAOH 

164 

Naa 

NA2 

165 

Na 2 0 

NA20 

166 

Ne 

NE 

167 

0 

0 

168 

OF 

OF 
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Species 

Chemical 

Key Punch 

Number 

Name 

Name 

169 

OF a 

OF2 

170 

OH 

OH 

171 

Oa 

02 

172 

P 

P 

173 

PCI 

PCL 

174 

PC1 3 

PCL3 

175 

PF 

PF 

176 

pf 3 

PF3 

177 

ph 3 

PH3 

178 

PN 

PN 

179 

PO 

PO 

180 

PS 

PS 

181 

P 2 

P2 

182 

P 4 

P4 

183 

S 

S 

184 

sf 3 

SF6 

185 

SO 

SO 

186 

so 2 

S02 

187 

S 2 

S2 

188 

S2F2 

S2F2 

189 

SI 

SI 

190 

SiCl 

SICL 

191 

SiCk 

SICL 4 

192 

SlF 

SIF 

193 

SiFCla 

SIFCL3 

194 

SiF Cl 

2 2 

SIF2CL2 

195 

SiF.Cl 

SIF3CL 

196 

SiF 4 

SIF4 

197 

SIH*. 

SIH4 

198 

SiO 

SIO 

199 

Si0 2 

SI02 

200 

SiS 

SIS 

°01 

2 

SI2 

202 

SnBr 2 

SNBR2 


SnCU 

SMCL4 

204 

UF t , 

UF6 

205 

Xo 

XE 

206 

Zn 

ZN 
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5.3.12 Subroutine C'DES 

This is the m*tn program for 0DE and corresponds to the Main Program de- 
scribed in Reference 1. Generally, the routine oerforms the following functions: 

1. Reads code cards THERM0, REACTANTS, 0MIT, INSERT, and NAMELISTS 
and directs flow of program accordingly. 

2. Stores THERM0 data on tape. 

3. Calls subroutine REACT to read and process REACTANTS cards. 

4. Reads 0MIT and INSERT cards and stores species names. 

5. Initializes variables in namelist $0DE. 

6. Reads and writes namelist $0DE. 

7. Converts assigned densities, if any, (RHO(i) in $0DE) to specific 
volumes: VLM(i) = l/RHO(i). 

8. Stores the number of pressures or volumes in NP. 


9. Stores values of o/f in 0XF array. If o/f values have not been input 
directly, they are calculated as follows: 


Values 

Code 

o/f calculation in 
main program 

Oxidant to fuel weicnt ratio, o/f 
Fuel to air weight ratio, f/a 
Percent fuel, %F 

Equivalence ratio, r 

Not specified 

OF = .TRUE. 

FA = .TRUE. 

FPCT = .TRUE. 

ERATI0 = .TRUE. 

0XF(i) = o/f 

0XF(l) = 1/ (f/ a) 

0XF(i) = (100-%F)/(%F) 

_ rV -(2) v +(2) 
rtvp/i) _ T v “V 

V^v 1 / _(1) + (1) 

rV U '+V u; 

rtvrfO - WP W 

0XF(l) WP(2) 


Values of WP(1) and V/P(2) are defined in appendix B or Reference 1. 

10. Makes necessary adjustments to consider charge balance if I0NS= 
.TRUE. . This is done by adding 1 to NLM and E to LLMT array. 

11. Calls SEARCH to pull required THERM0 data from tape and to store the 
data in core. 

12. Sets initial estimates for compositions. These estimates are set with 
each $0DE read. They are used only for the first point in the lists of 
variables in namelist (e.g., the first o/f and the first T and P in a TP 
problem). /II succeeding points use results from a previous point for 
estimates. 
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For the first ooint the program assigns an estimate of 0.1 for n, the total 
number of kilogram -moles per kilogram. The initial estimate of number of moles 
of each gaseous species per kilogram of mixture n^ is set equal to 0.1/m where 
m is the total number of gaseous species. Condensed species are assigned zero 
moles . 

13. Sets IUSEft) positive for condensed species listed on INSERT cards 
(see IUSF. array). 

14. Calls R0CKET is RKT is true. 

Note: while the use of oxidizer to fuel ratio, 0/F, is not generally used 
in describing solid rocket motor propellants, it can be used to great advantage in 
specifying binder, metal loading, and oxidizer ratios in parametric studies. 


5.3,13 Subroutine OMEGA 



! 


I 

’ 

» 

% 

l 

| 

> 


This subroutine calculates the exponent, used in the viscosity-temper- 
ature relationship 

T O) 

(r-) 


u (i ref v T 


ref 


using the method of least squares. That is, it calculates the value of fJt3 which 
gives the smallest sum of the errors squared. This form of the viscosity- temper- 
ature relationship was selected since both the TD2P and TBL modules require 
viscosity data in this manner. 

In order to supply the maximum amount of accuracy and also to minimize the 
variation in data due to the selection of an exit area ratio, it was decided to match 
the throat value of viscosity exactly and select an which would provide the best 
fit for viscosity at the chamber and exit of the motor. 

The form of the error, E, was taken to be 


E = In a/fi* " co l n T/T* 

Squaring the errors, differentiating with respect to oo , and setting the results 
equal to zero, yields the following value for ^ 

os = (In T c /T* In 4 / M * + In T/T* In n f / fl *)Aln T/T*) 2 + (In T/T*) 2 ) 


where 


T = temperature 
H = viscosity 
e = refeis to the 
c .= refers to the 
* = refers to the 


exit plane 
chamber 
throat plane 
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5.3.14 Subroutine 0UT1 

This subroutine, together with entries 0UT2 and 0UT3, writes statements 
common to all problems. 0UT1 writes statements giving the data on REACTANTS 
and on o/f, percent fuel, equivalence ratio, and density. 

Entry 0UT2 . - This entry writes the statements for printing values of pres- 
sure, temperature, density, enthalpy, entropy, molecular weight, fe In V/d In P),^ 
(if equilibrium), fe In V/fe In T)p (if equilibrium), heat capacity, yg, and sonic 
velocity. These variables and corresponding labels are printed with a variable 
format described in BL0CK DATA. 

Entry 0UT3 . - Entry 0UT3 writes statements giving the equilibrium mole 
fractions of reaction species. 

Subroutine 0 CTT1 has also been modified to print out the following addition- 
al items: 

a) gas velocity 

b) molecular weight of the gas-condensed phase combined 

c) mass fraction of condensibles 

d) molecular weight of the condensibles in the chamber 

e) molecular weight of the gas in the chamber 

f) the erosion parameter, 

In addition, 0UT1 also traps the starting conditions for the 0DK module. 
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5.3.15 Subroutine REACT 


The purpose of subroutine REACT Is to read and process the data on the 
REACTANTS cards. The subroutine is called from the main program after a REACTANTS 
code card has been read. The data on these cards are described in the REACTANTS 
Cards section (p. 62) of Reference 1 and in Volume III, Section 2.6.1. References 
to page numbers and equations given below also pertain to Reference 1 . 

The reactants may be divided into two groups according to card column 72 
on th? REACTANTS cards. The two groups are oxidants (O in cc 72) and fuels (cc 
72 £ O). We generally keypunch F in card column 72 for fuels even though this 
is not necessary. The contents of card column 72 are read into F0X. Depending 
on the contents of F0X, program variables relating to oxidants or fuels are sub- 
scripted 1 for oxidants and 2 for fuels. 

The FORTRAN symbols for the properties read from the REACTANTS cards 
and their associated properties (discussed in INPUT CALCULATIONS, p. 55 of Ref- 
erence 1) are as follows: 


Property 

FORTRAN symbol 

a (k) 

a ij 

wj k) 

Nj k) 

(H° ) (k) 

(u^) (k) 

j 

00 

p j 

ANLM(j,m) a 

PECWT(j) (if no M in cc 53) 

PECWT(j) (if M in cc 53) 

ENTH(j) (if not UV Problem and 00 not in cc 37 and 38) 

ENTH(j) (if UV problem and 00 not in cc 37 and 38) 

DF,NS(j) 

1 


a Each of the j REACTANTS cards contains from 1 to 5 stoichiometric coefficients 
read (indicated by subscript m) into ANUM(j,m) and their corresponding chemical 
symbols read into NAME(j,m). In relating an ANUM(j,m) with afj^ , the index i 
associated with a particular chemical element is determined from the chemical 
symbol in NAME(j,m). 
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If there are several oxidants their properties are combined by subroutine 
REACT into properties of a total oxidant using the relative proportion of each 
oxidant given on the REACTANTS cards* Similarly, if there are several fuels, 
their properties are combined into properties of a total fuel. The total oxidant 
and total fuel properties discussed in INPUT CALCULATIONS^ and their associated 
FORTRAN symbols are as follows: 


Property 

FORTRAN symbol 

Equation 
(see Ref, 1) 

b[ k > 

BCP(i,k) 

(187) 


Rr<fW(j) 

(190) 

v (k) 

HPP(k) (if not UV problem and 00 not in cc 37 and 38) 

(192) 


HPP(k) (if UV problem and 00 not in cc 37 and 38) 

(194) 


AM(k) 

(196) 

(k) 

0 

RH(k) 

(198) 

V + ( k) 

VPLS(k) 

(200) 

v” M 

VMIN(k) 

(201) 


PECWT(j) 



If any of the p are zero then RH(1) = RH(2) = 0. 

These total oxidant and total fuel properties are subsequently combined 
into total reactant properties by using the values of oxidant-fuel mixture ratios 
obtained from the main program. This Is done in NEW0F, an entry point in SAVE. 

Other common variables set by REACT are LLMT, NAME, ANUM, ENTH, FAZ, 
RTEMP, F0X, DENS, RMW, M0LES, NLM, NEWR, and NREAC. 

A provision is made for eliminating a second tape search when two consecu- 
tive sets of REACTANTS cards contain the same elements. This is done by saving 
the element symbols (LLMT( ?)) in LLMTS(i), the kilogram-atoms per kilogram 
(BOPU ,k)) in SB0PU,k), and the number of elements (NLM) in NLS, 

Atomic weights used in equation (190)^ are stored in AT0M(2,i). The 
corresponding chemical symbols are stored in AT0M(l,i). The oxidation states 
of the chemical elemenis V*or V~ used in equations (200) and (201)^ are stored 
in AT0M(3,i), The AT0M array is stored by BL0CK DATA. 
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5,3. 16 Subroutine RKT0UT 


This subroutine calculates various rocket performance parameters from pre- 
viously calculated thermodynamic properties. 

It is also the control program for writing rocket performance output. It con- 
tains the WRITE statements that apply specifically to rocket parameters and it 
calls subroutine 0TJT1 and entries 0UT2 and 0UT3 for the WRITE statements common 
to all problems. The rocket parameters are printed with the variable format, FMT, 
described in BL0CK DATA. 

Subroutine RKT0UT is called from subroutine R0CKET. 
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5.3.17 Subroutine R0CKET 

This subroutine is the control program for the RKT problem (rocket perfor- 
mance calculations discussed in section R0CKET PERFORMANCE).^ A flow 
diagram for this subroutine is given in Figure 5 of Reference 1 0 Subroutine 
R0CKET obtains the required thermodynamic properties for equilibrium performance 
by calling subroutine EQLBRM. For forzen performance, subroutine R0CKET calls 
subroutine FR0ZEN to obtain the required thermodynamic properties. Rocket per- 
formance parameters are then obtained by calling subroutine RKT0UT. In addition 
to calling RKT0UT, and FR0ZEN, and in addition to using controls common to all 
problems (discussed in section MODULAR FORM OF THE PROGRAM, p. 75, Ref- 
erence 1) subroutine R0CKET also does the following: 

1. It calculates estimates for throat pressure ratios. 

2. It calculates estimates for pressure ratios corresponding to assigned 
area ratios (if any). 

Subroutine R0CKET was modified to perform the following tasks: 

a) compute the area ratio at which solidification of condensed phases 
would occur if condensed phases are present. 

b) on option, restricts the total amount of condensed phases during 
the expansion to the amount Initially in the rocket motor chamber. 

c) computes the print out station at which subroutine OUT1 will trap 
the start conditions for a kinetics, 0DK, calculation. 

d) computes, using subroutiens MATCH, MUK, and 0MEGA, transport 
properties for the BAL, TD2P, and TBL modules. 

e) writes (punches) out the linkage data necessary for communications 
with the various loss modules. 
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5,3,18 Subroutine SAVE 

This subroutine has several functions, all of which are concerned with sav- 
ing some information from a completed calculation for subsequent use in later cal- 
culations, The primary purpose is to save computer time by having good initial 
estimates for compositions. 

These estimates for the next point, NPT, come from either the point just 
completed, ISV, or some other previous point. The flow of the routine is directed 
by ISV as follows: 

1. ISV positive. Transfer compositions from the point just completed for 
use as initial estimates for next point (transfer EN(j, ISV) to EN(j, NPT)), 

2. ISV negative. Save values of ENLN(j) for gases and EN(j) for condensed 
in SLN(j) , ENNin ENSAVE, ENNL in ENLSAV, IQ1 in IQSAVE, JS0L in 
JS0LS, JL1Q in LKIQS, and NLM in LL1. (These values are saved because 
they are to be used as Initial estimates for some future point and they may 
be overwritten in the meantime.) Make ISV positive and transfer EN(J,ISV) 
to EN(j,NPT). 

3. ISV zero. Use the data previously saved (as discussed in 2. as initial 
estimates for current point. Restore IUSE codes and inclusion or exclu- 
sion of "E" as an element for I0NS option. 

Entry NEW0F . - NEW0F combines the properties of total oxidant and total 
fuel calculated in subroutine REACT with an o/f value to give properties for the 
total reactant. NEW0F is called for each mixture assigned in the MIX array in 
$0DE namelist. It is called from subroutine R0CKET. The calculated properties 
and corresponding FORTRAN symbols are as follows: 


Property 

FORTRAN symbol 

b °i 

h(/R 

u^/R 

p o 

r 

B0(l) 

HSUBO (If not UV problem) 

HSUBO (If UV problem) 

RH0P 

EQRAT 




Equation 


(191) 

(193) 

(195) 

(199) 

(204) 


Subroutine HCALC Is called by Entry NEW0F to calculate the enthalpies 
for each reactant that has zeros keypunched In card columns 37 and 38 In its 
REACTANTS card. 

Values of HPP(2), HPP(l). HSUBO, B0P(l,2), B0P(t,l), and B0(l) are printed 
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5.3.19 Subroutine SEARCH 


This subroutine selects the Thermodynamic Data to be used in the problem. 

A scan is made of the master Thermodynamic Data tape and those species that are 
consistent with the chemical system under consideration are selected. As the 
thermodynamic data are being selected, the subroutine also complies a set of 
formula numbers, a^, from the formulas of the reaction products. A short Thermo- 
dynamic Data file is also generated for use in subsequent calculations, 

A check is made near the beginning of the routine to prevent THERM0 data 
from exceeding their storage allotments. These variables are all in labeled common 
SPECIES and are currently dimensioned for 150 species. However, this dimension 
may be reduced to save storage. 

SEARCH is called from the main program when the logical variable NEWR 
is true. NEWR is set true in REACT to indicate a new chemical system. REACT 
also stores chemical element symbols for the current chemical system in the LLMT 
array. SEARCH stores THERM 0 data in core for each species whose elements are 
included in the LLMT array (unless the species name was listed on an 0MIT card). 

The THERM0 data are stored in common variables TL0W, TMID, THIGH, 

SUB, A, C0EF, and TEMP. SEARCH writes out the names and dates of species 
whose data are stored in core. 

SEARCH initializes the IUSE array. IUSE(j) for gaseous species are set 
equal to zero. IU3E(j) for condensed species are set equal to negative integers. 

For the chemical system under consideration the first possible condensed species 
is set equal to -1, the second to -2, and so on, with one exception. In the event 
there are two or more condensed phases of the same species, set equal to -4, for 
example, IUSE(j) for B^O^s) will also be set equal to -4. A description of the IUSE 
array is given below. 

The various condensed phases of a species are expected to be adjacent in 
the THERM 0 data as they are read from tape. These phases must be either in 
increasing or decreasing order according to their temperature intervals. 

NS contains the total number of species stored in core. NC contains the 
total number of condensed species (counting each condensed phase of a species 
as a separate species). 
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IUSE array . - Each value in the IUSE array Is associated with a species. 
These values of IUSE serve two purposes: 

1. They indicate which species are to be included in the current iteration 
(IUSE (j) <0 for excluded species and IUSE(j) ^0 for included species). 

2. They indicate multiple phases of the same species if absolute values of 
IUSE(j) are equal. 

The IUSE(j) are initialized in subroutine SEARCH and the main program as 
follows: 

1. IUSE(J) = 0 for all gaseous species. 

2. IUSE(j) = n for all condensed species whose names have been listed on 
INSERT cards. The number n indicates the species was the n*h condensed 
species whose THERM 0 data were read from tape. 

3. IUSE(j) = ~n for all condensed species not listed on INSERT cards where 
n is defined in 2. 

These initial values of IUSE(j) may be adjusted later in subroutine EQLBRM. 
For condensed species, the sign is adjusted as species are included or excluded 
in the current iteration. 

For the I0NS option, IUSE(j) values for ionic species are set to -10000 when 

_o 

the mole fractions of all ionic species are less chan 10 




5.3,20 Subroutine SHCK 

Subroutine SHCK Is a dummy routine which replaces subroutine SHCK of 
Reference 1. 


5.3.21 Subroutine THERM P 

Subroutine THERMP Is a dummy routine which replaces subroutine THERM P 
of Reference 1. 


5.3.22 Subroutine VARFM T 

Subroutine VARFMT (variable format) adjusts the number of decimal places 
printed in F-format in the variable format, FMT, according to the size of the 
number. It Is used for P </? , P, and A e / A t - Variable format is described in 
BL0CK DATA , Section 5.1.2. 


5,4 Link 30 Subroutines 


The Link 30 subroutines control the grain design (GD) and ballistics (BAL) 
calculations made in the SPP code. These two modules have been integrated to- 
gether from the programs of Reference 2 and Reference 6. The remainder of this 
subsection gives an overall description of the integrated ballistics module. 
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5.4,0, 1 General 


The Interior Ballistics Module models the Internal ballistics of solid rocket 
motors. Including motor environmental effects on burning rate (erosion, radiation), 
effects of radial slots between segments on the port flow, port flow pressure drop, 
and delayed burning (ignition) in deep narrow slots. Semi-empirical correlations 
for combustion inefficiency and nozzle throat erosion also are included. This 
computer module was adapted from an existing model used to analyze 
nozzleless solid rocket motors (Ref. 6), because it has gasdynamic features 
lacking in many published internal ballistics programs. The intent was to predict 
internal ballistics for motors ranging from small tactical to large boosters wherein 
aspects of these features are encountered. The major effort in the present program 
development Involved updating certain model expressions, and rendering the com- 
putational procedure compatible with the Grain Design Module and the rest of the 
SPP code. The model will be subject to continued revision as new information be- 
comes available. 

The logic of solving the ballistics equations is presented in the following 
subsection. A discussion of the module subroutine structure is presented after that. 
Analogous discussion of the grain design subroutine structure is contained in Ref. 

(2) and need not be repeated here. Finally, the details of program Input and out- 
put, for the integrated subprogram, are given in Volume III, The Program User's 
Manual. 

The overall geometry-aerodynamics interface between the GD&BM is the 
generation of port area change as a function of length in the Grain Design routines 
and return of web burned from BM to GD module, at the same length stations. 
Ballistic calculations may be performed for one or more time increments within an 
increment of bum (NB) specified in the GD input, but as the increment is predicted 
to bum out the actual web expended is returned as the starting point for the next 
increment of geometry calculation. 
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5«4»0.2 Solution Logic 


The general scheme employed in solving the equations of one-dimensional 
flow is shown in simplified terms in Figure 5.4-1. The procedure that follows 
the initialization of a case is as follows: 

• A head-end pressure is assumed. The calculations for a single axial in- 
crement are made and iterated to completion. The successive increments 
down the port are then calculated until a match point on the nozzle entrance 
is reached. * 

• At the match point, the local static pressure is checked against the pres- 
sure from throat flow constraints. If the two values are not within tolerance, 
a correction to head-end pressure is made and the calculation is restarted 

at the head. As successive iterations are made to satisfy the choke con- 
straint (outer loop), the axial stations (inner loop) are continuously iterated. 

• Parameters used in transient terms and integrating calculations are saved, 
and the time increment is advanced. The calculations resume with a new 
head-end pressure assumption. 

• As the calculations reach the end of a specified run time, or the motor 
has burned out, output is furnished to the Nozzle Performance Module. 

The computational logic, thus, is arranged as a triple-loop structure with 
the axial increment (inner loop), choke constraint (outer loop) and time increments 
(stepping). 

A more detailed diagram of program logic is shown in Figure 5.4-2. This 
representation is still simplified with regard to the individual branching for special 
cases, but is intended to describe the function of groupings of program steps in- 
ternal to the program, and presumes valid input. Figure 5.4-2 is organized in the 
same order as the source program. 

Ihe compiled-in constants are established and some core areas are cleared 


to ze r o. 


A subroutine to control the reading and organizing of tabular data is called 
b\ Subroutine READIN on completion of geometry inputs. The remaining data (in- 
tegers and real data ir.out) are returned in intermediate storage which is inspected 
for non-zero (nz) values. Any nz values are then placed in working storage for 
use hi MAIN 3. 


For an end-burner, the ‘head end" will be the grain face. A protective statement 
was included providing that if the motor port/throat is less than 1/4, the calcula- 
tion will proceed to the next increment. Thus, the end face is eventually located. 
Internal burners are unlikely to have such low port/throat, so this provision should 
not introduce a problem in internal burners. 
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Hie outside diameter or periferal geometry of the grain is established from 
tabular inputs which are fed to subroutine L00PS in the GD program. The GD pro- 
gram initializes the grain and computer bumback geometires. 

The following groups of calculations are repeated at each time increment 
of the problem. The igniter flow, if any is determined, and time dependent values 
of Aj# Ae, thrust, ambient pressure and thrust coefficient are determined either 
from tabular input or internal feedback. Calculations are made for each axial 
increment in the motor. The axial increments correspond to the GD X-grid 
plus an aft plenum, plus the throat. Trial values for the initial time step are used 
from the previous axial station; all other trial values result from the previous iter- 
ation of the station. Special case logic for ignition options, or initial time are 
not shown in the figure or discussed here. (However, they should be recognized 
in the program listing as IGN.GT.O, or HME.EQ.O. branches.) 

The port geometry is calculated by the GD module from the web burned. 
Results are interpolated within the bum interval to determine local bum volumes. 

• Base strand bum rate for the propellant at the axial station, as a function 
of local pressure, is found from an input table. The input data can be 
keyed to use log -log interpolation. The motor environment contributions 
are then calculated. 7 7 ^ effects are Incorporated by changing the base table. 

• If appropriate, logic and branching related to radial slots is performed. 

The logic path returns several places depending on the interpretation of 
radial slot flow. 

•Local conditions are calculated using the energy, state, and continuity 
equations, and thermochemical data. 

• Traps are executed to prevent crossing Mach 1, and the momentum equation 
is used to solve for pressure. The set of equations at each axial station 

is iterated with logic for subsonic flow or Mach 1, with generation of an 
error term to later correct head pressure; excess iterations will be shown 
as output. 

• Looping logic allows 10 iterations and then accepts the result whether 
converged or not. 

• The node number of the axial increment is checked to branch to the next 
node or to throat calculations. 

• Throat constraints are calculated from mass flow rate, nozzle throat area, 

C*, etc. 

•Nozzle entrance static pressure is calculated from throat stagnation pressure, 
assuming an isentropic expansion, to compare with the pressure determined 
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by Integration down the port. An updated head pressure Is calculated from 
the error In nozzle Inlet static pressure, or a term Is created to Indicate 
an overchoked prior calculation. 

*** i^ iroat constra i nt is not satisfied, loops are counted to 10. From 10 
to 20 loops, the enror Is printed; after 20 loops, the result Is accepted. 

The printed error messages Indicate whether the calculations are oscillating 
or are overdamped and allow the user to determine whether the run Is acceptable. 

• Values used In transient terms or Integrations are stored for use In the next 
time step. Parameters along the motor are printed. The burn-back of radial 
slots and wet are updated. Axial nodes are coded out Just prior to bum- 
back of the radial slot face to the node position. For radial slot combus- 
tion the slot ends are Input to the GD program as non-burning. The 

bum-back Is then explicitly calculated In MAIN3 and the axial Increment 
locations are modified. 

• Nozzle exit parameters and thrust are calculated from isentropic expansion 
equations, and a series of running Integrals are calculated and output. 

• The values of time and Input run time are checked for completion. The 
next time Increment or rewind follows this branch. 

The Indexing scheme employed In the program Is shown In Table 5.4-1 
as an aid In following the details of the logic. This logic Is designed to treat non- 
perpendicular bum-back of radial slots; but Is also used to "key" pressure drop 

calculations such as a sudden expansion or non-dlffuslng taper (l.e., separated 
flow). 


Table 5.4-1 INDEX SCHEME 


Symbol 


Description 


I 

Axial node location of calculation. 

I 

Prior axial node, usually J=I-1, J>0, but 
is varied to skip cast-off nodes (for example 
J=I-4 if three nodes are burned out). 

L 

Radial slot number, 1 at head end; slot re- 
lates to following segment. 

N 

Node number of nczzle entrance (36 in 
sketch) . 

NN 

Throat node. 


8 









Thrust calculations and output 
Aero calculations and output 
Ballistics calculations and output 
Integrals and output 


Check for completion 
of burn Increment 


I 

Check for completion 

No - Advance time Increment 
yes 

No 


yes 



RETURN KBAL 


RETURN TO GD 


Figure 5.4-2 Module Structure 
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Figure 5.4-3 Module Logical Structure 




Figure 5.4-4 Input Groupings Following Grain Geometry Cards 



5 ,4 . 1 Program LINK40 


This subprogram is mainly intended to facilitate the conversion of the SPP 
code overlay structure to other computer systems. Its only function is to call sub- 
routine BALIS. 

5.4.2 Subroutine BALIS 

This subroutine reads and writes the namelist data set $BAL. It then calls 
subroutine READIN and FIG1 to read in and then calculate the grain design parameters. 
Control of the ballistics and grain design calculations sequences are then handled 
in subroutines MAIN3, PVS and L00PS. MAIN3 performs the axial looping for bal- 
listics calculations, while L00PS does the same for the grain geometry calculations. 


5.4.3 BL0CK DATA 

This routine is used to zero out the common blocks used in this module. 
The KALPHA array in labeled common C0DE is also initialized. 

5.4.4 Subroutine C0NTR1 

This subroutine has been written for use in the Grain Design Program as 
an auxiliary program for CONTR2. The routine computes the initial data needed 
by CONTR2 that is dependent only on the type of cone (or cylinder) being consid- 
ered and independent of the burn distance, X plane, or line in question, 

5.4.5 Subroutine C0NTR2 

This subroutine has been written for use in the Grain Design Program. 

The purpose of the routine is to determine the points of intersection, if any, 
made by a line and a frustrum of a cone (or cylinder) located in three dimensional 
space. 

5.4.6 Subroutine C3TEFF 

This subroutine calculates the empirical c* efficiency described in Volume 
I, Table 4-2. 
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5,4.7 Subroutine DA TAIN 


TWAIN is n random rnnd routine tint Interprets input of selected variables. 
UATA1N operates under the control of TABIN, which requests the reading of cards 
with integers or real data, which is in turn called from Subroutine READIN. 

5,4.8 Subroutine ERROR 

This subroutine was written for use with the Grain Design Program. The 
routine checks the input data for errors and writes error messages. See Appendix 
A for conversion aides. 


5.4.9 Subroutine FIG1 

Subroutine FIG1 has been written for use with the Grain Design Program. 
The program calls the figure routines for initial calculations. The initial bum 
used to input geometric figures with rounded comers is monitored by subroutine 
FIG1. 


As Subroutine FIG1 is entered, I is initialized to one for the first figure. 

ID is computed from IDENT(I) and a "computed go to" is executed on ID in order 
to call the proper figure routine (CONTRJ or PRISMl). On return from the 
figure routine, I is incremented by the proper amount and a test of I against 
NOSRF is made to determine if more figures remain. Upon completion of the figures, 
control is returned to the main program. 

5,4.10 Subroutine FIG2 

Subroutine FIG2 has been written for use with the Grain Design Program. 

The program calls the proper figure routines for computing the points defining the 
void volume of the propellant. 

As Subroutine FIG2 is entered, the motor case boundaries are computed for 
the specified Y value being considered. Then the bum distance to be considered 
(TRIAL) and the identification code (ID) for the figure to be considered are com- 
puted. Then the proper figure routine is called and upon return Subroutine UNIOr 
is called to form a union between the computed Z points. These points represent 
the shape of the propellant of the motor case. 
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If all input figures have not beer covered, the loop is repeated, from the computing 
of TRIAL, for the next input figure. 

5.4.11 Subroutine HNDX 

This routine performs linear interpolation or extrapolation in tables while 
saving its place in the table of independent variables. 

5.4.12 Subroutine L0(2flPS 

Subroutine L00PS sequences the grain cfesicn subroutines through the input 
lists of Y increments and performs a summation cf port area contributions for the 
Ntr axial station in the motor. This subroutine replaces several in the original 
program of Pef. (2) since less options are retained for ballistics calculations and 

inverted order of sequencing is not needed. Thus, in the revised version, the outer 

? 

loop is always bum increment; the middle loop is always motor axial station (di- 
mension x); and the inner loop is on the radial intersection (dimension y). The sub- 
routines performing the actual geometry calculation are essentially unchanged from 
Ref. (2): changes were introduced in the sequencing routines to imporve program 
run times . 

5.4.13 Subroutine MAIN3 

This routine controls the calculation of both the ballistics and grain design 
modules. In fact it is the main portion of the ballistics calculation. 

Subroutine MAIN3 replaces a series of subroutines controlled by MAIN3 of 
the original Ref. (2) program. This subroutine calculates internal motor ballistics 
based on grain geometry output from the grain design subroutines. Control is 
retained in MAIN3 fcr as many time increments of web burned as necessary to pre- 
dict bumcut of the particular increment of web calculated by the grain design pro- 
gram. The web calculated by the ballistics subroutines is then returned to the grain 
design program through subroutine BALIS for use as the starting web in the next 
Increment. 

An artificial point in the nozzle entrance region is introduced as a match 
point for calculations proceeding down the port to the nozzle, for the purpose of 
testing the choking constraint without having to calculate flow near Mach one. 
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The area of the "nozzle entrance" is not critical to program operation, but run time 
is reduced if a large area is selected. The pressure computed by proceeding down 
the port must agree with that computed (by assuming isentroplc f;ow) from the 
throat back to the artificial point, within a certain tolerance. 

Values of c* efficiency and throat erosion, which are output from the bal- 
listics subroutines, are based on semi-empirical correlations. Thrust may be cal- 
culated by the ballistics subroutines standing alone, but requires an input nozzle 
efficiency. This nozzle efficiency may be calculated external to the computer por- 
gram, using such empirical expressions as contained herein or elsewhere. A 
series of motor operating parameters and running integrals are output and saved 
for use in other modules. If the other modules are exercised, the thrust parameters 
calculated herein are overridden. 

Detailed output of conditions at each axial node in the calculations may be 
requested as a function of time using tabular input of the increments desired. 

Calculation of the bum-back of radial slots is made by: 

1) entering the slot faces as "nonbuming" in the grain design program 
input; 

2) overriding the grain inner diameter calculation with the outer dia- 
meter in the slot; and 

3) explicitly calculating the location of the end of the port in sub- 
routine MAIN3 

The node at the end of the grain retains its identity number but is then a moving 
node. As the burning face approaches the next x-node along the grain, the in- 
terior node is codec out. Table input is provided to allow for ignition delay in the 
slot (radius ignited versus time). If significant delay is used, the face of the grain 
will not remain perpendicular to the motor axis (axial position is retained at the 
inner diameter); additional table input is then needed to correct the grain length 
as the port bums outward. Evidence that corrections are needed is that the total 
propellant weight expended will be under-predicted. The burning rate on the slot 
face is assumed to be the input strand rate. 

5.4.14 Subroutine NER0DE 

This routine calculates the nozzle throat erosion rate from input tables 
or from a correlation determined from analytical calculations. 
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5,4,15 Subroutine PRISM 1 


This routine calculates the initial data needed for each specific prism. 

Only right triangular prisms are considered. 

5,4,16 Subroutine PRISM2 

This subroutine has been written for use with the Grain Design Program. 

The purpose of the routine is to determine the two points of intersection, if any, 
made by a vertical line and a triangular prism located in three-dimensional space. 

When a prism bums out, the comers bum into spheres, the edges bum 
into cylinders, and the end and side planes bum straight out. A line will inter- 
sect a prism either through the cylinder, sections (edges), sphere portions (cor- 
ners), or not at all. 

This routine considers the intersecting line as always a vertical line in 
the standard three-dimensional system. Hence, the vertical line is defined by 
an X and a Y vlaue. The number of line intersections necessary depends on the 
particular prism beirg considered and the accuracy desired in subsequent computa- 
tions. 

Input to this subroutine consists of an X value, a Y value, a T (bum) dis- 
tance (positive for out-burning, negative for in-burning), and those computations 
made in PRISM1. 


5.4.17 Subroutine PVS 

This subroutine is used only with the ballistics mode. It provides peri- 
pheries, volumes, flow areas, and burning surfaces areas (and totals) at or be- 
tween each X-station and bum interval (time, to time .). 

K rC T 1 

Cross sectional flow areas at time k and time k+1 are supplied to this ro- 
utine for each X-stat*on. These are provided in two tables thus: 

Aj j = Cross sectional flow area at time k and X^ . 

A 2 j = Cross sectional flow area at time k+1 and X^. 

The n X-stations' locations are provided in a table of X^. 
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Volume burned between each two adjacent X-stations during time, to 

Jv 

time k+l com P utec ^ and summed: 

V total = £ AVj 


A V j f (A 2,j " A l,j* + ^2,j+l * A l,j+1^ ' ,X j+l _X j! 

The average periphery of the flow area between time^ and time^ + j , at 
Xj is calculated as the change in flow area over distance burned at that X-station. 
'Diis is calculated only to be printed. 


P = _U. 1J. 

J ABj 

The average burning surface area during the interval time^ to time k+ j, 
and X^ to X^ +1 is calculated and summed,, For this calculation, the equation 
S=^V/^B had to be modified because is not necessarily normal to S in the 
ballistics mode. 

5.4.18 Subroutine READIN 

Subroutine READIN has been written for use with the Grain Design Program. 
The program initializes variables and reads the data. A test is made on the data 
in order to pass control to the appropriate place for storage of the data. 

5.4.19 Subroutine SFERE2 

This subroutine has been written for use with the Grain Design Program. 
The routine will determine the intersection of a line with a sphere located in 
three-dimensional space. 

The routine solves the equation 

C = R3 - (Y c - - (X c - X 2 )a 

where R = The radius of the sphere. 

X^ = Denotes the location of a plane in three-dimensional 
space. 
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Y^, = Denotes a line in the X plane. 

Y = The Y coordinate describing the center of the sphere. 
Xj = The X coordinate describing the center of the sphere. 
If C is negative, the line does not intersect with the sphere. 

If C is zero, the line is tangent to the sphere. 

If c is positive, the intersection points of the line with the sphere are given by: 

Z = z x + C 

where \ = The Z coordinate describing the center of the sphere. 

5.4.20 Subroutine STORE 

Subroutine STORE has been written for use with the Grain Design Program. 
Its purpose is to determine the nature of the input data and store it in locations 
to which the main program has access. 

Control is sent to the program from Subroutine READIN by the statement 
CALL STORE (J,JO); JO has been assigned an integer value of 1 to 4 by Subroutine 
READIN depending on the card type and is used in a "computed go to". J is a 
subscript incremented in Subroutine READIN for the correct storage of data into an 
array. 

5.4.21 Subroutine TAB 

TAB is a routine to interpolate (or look up) values from the tabular inputs. 
Linear interpolation or log-log interpolation is used and no extrapolation is made; 
rather, the value at the boundary of the table is returned. The tables are assumed 
of the form 

z = f (x,y) 

with a call statement of the form 

CALL TAB (X, Y, Z, No) 

where No is a table number assigned in Subroutine MAIN3. The size, shape, type 
of interpolation, and order of table loading are established by input of the tables 
and thus are not transmitted to other subroutines. 

Internally, the logic is established to search the table by stepping from the 
prior point to the new one, consistent with an iterative program with highly re- 
petitious call statements, in contrast to computing indexes. 
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5,4.22 Subroutine TABIN 

TABIN is a routine which prepares input of tabular data. The addressing for 
the tables is arranged, indexes are initialized, and limits to table size established. 
Input of integers and non-tabular program inputs are simply passed back to READIN 
for interpretation. 


5.4.23 Subroutine UNI0N 

This subroutine has been written for specific use with the Grain Design 
Program. The subroutine performs a union with Z points solved for previously and 
those Z points already entered in the table. 


A table of void points or grain points is kept by entering, deleting, or in 
general, performaing a union with paris of Z points. These Z points represent the 
intersection of a line and the basic figures allowed for in the Grain Design Pro- 
gram, Based on a flag, the Z points can be treated either as cuts or fills to the 
table. Only those points equal to, or contained within, the case boundaries are 
considered. All points in the table are left adjusted. A count of the number of 
pairs is also computed. 


Program Symbols 
Fortran 


Description 


TABLE 


NOZE 

RDIN(2) 

RDIN(3) 

Z(l) and Z(2) 
ZEE and ZEET 
ZF 

J 

I 

I UFA 

KALI 

IDIF 

KUFA 


Location where void or grain points are stored after 
union is made. 

Number or paris of Z in table. 

Lower boundary or case limit points. 

Upper boundary or case limit points. 

The void or grain points. 

Z(l) and Z (2) . 

Temporary stroage for one of void or grain points. 
Number of points in table. 

A count of number of points already looked at in table. 

Flag used to enter or delete Z points. This flag is 
also set outside of union to tell if table contains 
grains or void points. 

A counter used to aid in entering or deleting points 
from the table. 

A difference value used to left-adjust points in the 
table. 

A branch flag used within the program. 
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5.5 Link 40 Subroutines 


Overlay 40 contains the 0DK module which calculates the loss due to finite 
rate gas phase kinetics as incorporated into the methodology programmed in the 
SPP code. 

This overlay includes several sub-overlays which are executed in sequence, 
i.e. # Link 42 is executed after Link 41, etc. The subroutine descriptions included 
in this section are in alphabetical order within the individual overlay structure. 

All of the subroutines described in this section are the results of the modi- 
fication of the 0DK program described in Reference 3. While only certain of the 
subprograms of that reference have been modified, the entire pertinent body (i.e., 
0DK documentation) of Reference 3 is repeated here. 


5.5.1 Sub r outine LINK4 0 

This subroutine consists only of a call to subroutine 0DK and is used only 
to facilitate conversion of the program overlay structure to other computers. 


5.5.2 Subroutine 0D K 

This subroutine acts as the driver for the one dimensional kinetic expansion 
calculation (0DK). It calls in over A ays LINK41 through LINK43. Subroutine 0DK 
also writes (punches) out the linkage data dealing with the restricted equilibrium 
option of the 0DE module. 


5.5.3 Subroutine S TF 


This subroutine evaluates the thermodynamic functions Cp^/R, H^/RT, 
S^/R, fror curve fit coefficients. The subroutine uses the same procedure as 
subroutine CPHS. The additional functions dfCp^/dT and free energy, G^/RT, 
are also computed. The calculated functions are then converted to the internal 
computational units for use by the kinetic expansion calculations. Also, low 
temperature thermod/namic data from tables is used when required if the LTCPHS 
directive was specified in the input to the SPP code. 


5.5.4 Subroutine STOICC 

For each reaction this subroutine constructs two vectors of stoichiometric 
coefficients, one for reactants and one for products. Up to 10 reactants and 10 
products may be considered for each reaction. The total number of entries in the 
resultant linear reaccion table is 600, i.e., the sum of all stoichiometric coef- 
ficients cannot exceed 600. 


5.5.5 Program Link 41 

This subprogram is mainly intended to facilitate the conversion of the 
SPP code overlay structure to other computer systems. Its only function is to 
call subroutine 0DKINP which handles the input to the 0DK module. 
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5,5/6 Subroutine ECNV 


This subroutine translates a BCD string of characters into one 
floating point numeric value. E,I, and F formats are permitted with the 
result always a floating point number. It is called by subroutine REAXIN to 
decode numeric fields in the species and reactions cards. The subroutine 
is coded entirely in FORTRAN. A BCD string of blanks will result in a floating 
point zero returned value. 


5.5.7 Subroutine NUMBR 

This subroutine converts a one character BCD number to a FORTRAN integer 
number. The subroutine is coded entirely in F0RTRAN. 


5-62 


5.5.8 Subroutine CfoKINP 


This subroutine provides the input processing for the kinetic expan- 
sion calculation, ft performs the following functions: 


1) Variable Initialization to nominal values 

2) Calls subroutine REAXIN to input the reactions cards and species 
cards if necessary 

3) For an 0DE-0DK problem, calls subroutine SELECT to select 
those species to be considered for the kinetic expansion 
calculation 

4) Reads $0DK namelist input data 

5) Converts nozzle geometric parameters from input units: inches, 
degrees; to internal computational units: feet, radians 

6) Computes nozzle tangent coordinates using: 


r t = 1 + R d (1 - cos 6) 
x t = R d sin 0 


7) For conical nozzles, computes the axial coordinate for the exit 
station from the following relation: 


x exit 

x exit 



r t + x t • tan 6 
tan 0 


- ( 1 + R d -VEP) 


2 



x exit 

x exit 




<X A 


8) ' For conical nozzles, the internal axial print stations 
are computed using: 


X, = 


VaRPRNT(J) - r t + x t • tan 6 
tan 0 


y j* x t 


X, = 


R 2 - ’ 

d 


1 + R d - (ARPRNT(J) ) 


4 


Xj <X t 
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9) The sum of input or selected species concentrations is checked 
for unity ( + XMFTST, where XMFTST is an input number), and 
then normalized. 

10) 11 the input parameter RZN0RM is input, the input contoured 
nozzle table is normalized by RZN0RM. 

In addition to the above, the following modifications have been made to 
the 0DKINP routine. 

a) If the nozzle geometry has been input via the $GE0M namelist, 0DKENP 
traps and stores the variables from this input into the correct locations 
and performs the appropriate unit conversions. 

b) Was modified to accept the parabolic and circular arc options (IWALL=2 
and 3) from the $GE0M namelist. This was done by computing the 
nozzle wall points and derivatives at 20 equally spaced points and 
substituting them into the tables used in the contoured wall input 
option. 

c) Was modified to accept as input the quantities needed to calculate 
the zero particle lag flow equations. 
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5,5,9 Subroutine REAXIN 


This subroutine processes SPECIES, REACTIONS, and THIRD B0DY REAX 
RATE RATT0S input cards. Reference may be made to Volume III, Section 2.3.2, 
the Program User's Manual, for a complete description of input requirements. A 
table of all species appearing in the input reaction set is generated for further 
processing by subroutine SELECT if required. 


5,5,10 Subroutine SELECT 


This subroutine provides the interface logic required to select the 
minimum species list required for the kinetic expansion calculations. The 
subroutine is only used for the 0DE-0DX. interface. The list of all species 
appearing in the input reaction set is matched against the list of species con- 
sidered for the equilibrium calculation. All species which appear in both a 
reaction and the equilibrium calculation list are selected for the kinetic expan- 
sion calculation. If a species appears in the reaction set but has not been 
considered for the equilibrium calculation, the program prints an error message 
and terminates the current case. 

If the INERTS directive was specified those species specified under 
that directive will be added to the list for the kinetic expansion calculation along 
with the other species selected. 

If the INERTS directive was not specified, all those species, con- 
sidered for the equilibrium calculation, whose mole fractions are greater than 
or equal to an input selection criterion will also be selected for the kinetic 
expansion calculation. Species selected in this way will be listed as inert 
species on the program output since they do not enter into chemical reaction. 

This routine was modified to select pairs of condensed phases if 
either of the phases passes one of the above selection criteria. 
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This subprogram is mainly intended to facilitate the conversion of the SPP 
code overlay structure to other computer systems. Its only function is to call 
subroutine PACK. See Section 5.5.13 below. 
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5.5.12 Subroutine CQfNVRT 


This subroutine converts input data from the externally input units 
to internally used computation units. In order to conserve computation time 
during the kinetic expansion, parameters such as molecular weights, are in- 
cluded in these conversions. Primed numbers are input quantities. 


a) Reaction rate ratio input for reactions requiring third body terms 
units: unities s 

internal units: ( lb s-mas s/lb -mole)" 1 
formula: XMM. . = XMM. ./Mw 

*■ 

b) Pre-exponential reaction rate parameter 
input units: cm, °K, g-mole, sec 

Q 

internal units: ft , 0 R, lb-mole, sec 


a j = 


aj * (.0160183)** 1.8 n j 
n 

n 

i=i 


n M V /ij 


Where X depends on the order of the reaction. 


and .0160183 = 


3.531 • 10' 5 ft 3 

, 3 

1 cm 


1 g-mass 


2.2* 10 3 lbs-mass 


c) Exponential Term: 

input units: kcaJ/mole 
internal units: °R 
formula: b. =bj • 905.770 


«*iiere 905.770 


1000 cal 
1 kcal 


1 I 8°R 

1.98726 cal/mole- : K "OFR 
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<9 Equilibrium Constant Multiplicative Factor: 
Input units: not input 

o 

internal units: (Ibs-mass) - °R/ft 
formula: 

n 

i-i Mw . Vii 

DkTEF(J) = —i-i 

n 

II 0.73034 


where .73034 = 49,721.011- A - goun . d \ ^ 

(lbs-mole) R 


1 atmos 2 

68,059.59 poundals/ft^ 


e) Pressure: 

input units: PSIA 

9 

internal units; poundals/ft 
formula: P = P’ - 4633.056 
where 

2 

4633.056 = - 14 - 4 * n • 32. 174 -&— 0 
1 ft z sec 


f) 


The initia 1 reference enthalpy is computed using 



i=l 


Subroutine C0NVRT has been modified to calculate the total amount of con • 
densed phase present at the kinetic expansion initial conditions. 


5-68 


$>5,13 Subroutine PACK 


On the basis of those species currently being considered, this sub- 
routine packs species and reaction information from the master tables into 
those control sections utilized by the one-dimensional kinetic expansion sub- 
program. 

The following is a sequential description of the packing procedures: 

1) Thermodynamic data for the species being considered is 
reed into core storage. 

2) The chemical species'molecular weights are computed 

3) The symbolic reactions are checked for mass balance. 

4) For a contoured nozzle the slope at each input wall point is 
computed using subroutine SLP. The wall coordinates, and 
each computed slope are printed for each input wall point and 
the print stations are set to the input axial coordinates. 

In addition to the above, the PACK routine processes condensed 
phases by setting R^O for those species and storing pointers up to 10 condensed 
phase species (5 pairs). This routine has also been modified to calculate the print 
positions for specified area ratios when the contoured (spline) wall option is 
selected. 
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5.5.14 Subroutine PRES 



This subroutine is used (when JPFLAG s 1) to compute the derivatives 
of an input pressure table. 

This subroutine is also used (when JPFLAG = 0) to generate a pressure 
table through use of an average expansion coefficient, Ne. The generated table 
extends from the initial contraction ratio through the nozzle attachment point plus 
one normalized throat radius. 

Input Pressure Table Derivative Computation (TPFLAG = 1) 


If a pressure table of NTB entries is input, the table of first deriva- 
tives is computed using: 



dP 

dx 



dP 

dx 


X NTB 


0 


P 


P 


(x n-H } ~ P(x n-1 } 
x (n+l)~ X (n-1) 

<»„> : Vi > 

x (n) " x (r«-l) 


, 1 < n < NTB 


, n = NTB 


The pressure at the initial axial position is obtained by interpolation 
using subroutine SPLN. 
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Internally Computed Pressure Table Computation (IPFLAG = 0) 


An average equilibrium pressure expansion coefficient from the 
chamber to the throat, is commited by iteration using subroutine 
SUBNE. The initial value for N (1 ' is 1.2. 


The approximate equilibrium contraction ratio at the initial axial 
position is computed from: 




where 


= pressure at the initial axial position 

P = equilibrium chamber pressure 
c 


A check is then made to determine the compatibility between the 
nozzle geometry and the requested contraction ratio. 


ifcT* 1 + [ R u +R iH i -cos %]• 

the circular arcs R u and R^ overlap and the following error message is printed: 

INLET GEOMETRY INCOMPATIBLE WITH INITIAL CONDITIONS 
The program will proceed to the next case. 
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Tables for pressure and its derivatives are constructed as functions of area 

ratio, a, and expansion coefficient, Ne. Formula 

used for the j + l iteration for pressure is: • * 


V 1 


P 0+1) p(l> 


fi-rfi ■ 


-1 


N -1 
e 


[2/(N +1)] 


(N e +l)/CN e -l) 


2/N 
(J)' 6 


1 - 




<N e -l)/N e 


- 1/2 


a -i 


for J = 1 


,( 1 ) : 


n+l 


p_! d(P/P c ) 

P c I + 

'x 

n 


(x 


n+l 


V 


where n refers to the n ^ table entry. 


The pressure derivative formula used is: 


d(F/P c ) 

dx 



5-72 


an 
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.r 


Next, tables for pressure and its derivatives are constructed 
by the program. Table entries are at increments of 


-x/7S 

for 

Xj < x < 0 

(R d sin e)/25 

for 

0 <x < R^ sin 0 


and 


1/25 


for 


R^ sin 0 < x < sin p + 1 


where the initial nozzle axial position, x 1# is computed from: 


- 

•\fei~- 1 - (R + R=) • (1 - cos e.) 

c i = " (R u + R i ) ‘ sln 6 i + teTeJ 


See Figure 5.5-1 below: 



Area ratio and its derivative and (a and are found by the five formulae 

below: 


1) x < Xj + Rj sin 

js. . ~ 2 lx ~** 1 .vr 

«* 


2) Xj + Rj sin 0j < x < - R u sin Qj 


a = “ R i^ _ cos “ ^ x_x i “ R i sin e i^ tan 6 i] 

~3x” = “ 2 *>/«"• tan0 i 


3) -R u • sin 0 i < x < 0 


- [■ * r » (> -f- -%) 


[ R U - x j 1/? - 


lx r— 

Y' V a 

“ I l /rt " 


4) 0 < x < * sin 0 




da _ 2x <JT 
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5) R d • sin 0 < x ^ R d • sin e + 1 
for cone 

a * £r t + (x-x^ * tan 0 

jSS-s 2 • a • tan 0 
dx 



for contour 
a ° Y 2 


da _ , v dY 
*T " 2 ' Y 'aT 


Three special points are included in the pressure table. These are 
a point at x = Xj such that 

P . _ P JL 


d( p / p c) 

— U 


and two points at x = 0 such that 

p = I p* 

p 

c / equilibrium 


He_ . f. 


3N e" 1 
2(N - l) 


dx -fi* [ N e 


2 

+ 1 


with R* = R u and R* = R d , respectively. 
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The following Items are Input directly to the computer program as described 
in Volume III, Section 2.8.5, Figure 2-2, left to right. 


u 


RI 

RWTU 

THETAI 

THETA 

RWTD 
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The purpose of this subroutine is to supply derivatives for a tabulated 
function. The end point derivatives may be specified or are calculated inter- 
nally by parabolic interpolation. Interior point derivatives may be found by 
a cubic spline fit procedure. 

Calling Sequence: 

X is a table of Independent variables, x 

Y is a table of the dependent variables , y^ 

N is the number of entries in each of the tables X, Y, and YP. 
1= 1, ...N 

MFLAG this entry is a flag, m, such that 

m > 0 Implies x is equally spaced 

m < 0 implies x is not equally spaced 

|m| = 1 y' will be continuous 
|m| = 2 y* arid y M will be continuous 
YP is a table of the derivative, yj 
W1 working storage of length N 
W2 working storage of length N 


below: 


W3 working storage of length N 


IFLAG 


this entry is a flag, i, such that 

i = o implies value for YP(1) and YP(N) will be calculated 
internally by parabolic differencing 

i = 1 implies values for YP(1) and YP(N) will be input 


i= 1 


Method 


The cubic spline fit procedure utilizes the interpolation formula given 


3 2 

y - A(x - x ) + B(x - x ) + C(x - x ) + D 
7 ' o £ o o 

y* = 3A(x - x q ) + 2B(x - x q ) + C 
y n = 6A(x - x q ) + 2B 

The piecewise cubic fit to a tabular function by the above relations will yield 
a discontinuity in the second derivative y M , between adjacent fits of: 


1 




tion equal to zero so that the second derivative is continuous across Juncture 
points. As applied to a tabular function, the above procedure results in a set 
of linear simultaneous equations (tri-diagonal) to be solved for the y| , provided 
that values for y' at the end points are known. 
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5,5,16 Subroutine SUBNE 

Calculates the average equilibrium pressure expansion coefficient 

from the chamber to the throat by iteration from the following formula 
(Newton's method): 



where -1.2. 

P* is the equilibrium throat pressure 
P c is the equilibrium chamber pressure 

This. subroutine Is used by subroutine PRES. 
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5,5.17 Program Link 43 


This subprogram is mainly intended to facilitate the conversion of the SPP 
code overlay structure to other computer systems. Its only function is to call 
subroutine MAIN ID which handles the control of the integration of equations de- 
scribing one dimensional reacting gas fl ow thro ugh a rocket nozzle. 


5.5.18 Subroutine EF 

This subroutine computes equilibrium, constants, K 


j 


« 


Kj EK(J) 
also computed are 


DATEP(I) 
'p A 



dK. 


l 

dT 

DKT(J) 

where: 

Ft i 


Ht i 


DATEF (J) 





l ‘U + & R, 


n 

E 


Ht 


i 





= species free energy at the current temperature 
= species enthalpy at the current temperature 
= is defined in subroutine C0NVRT, Section 5.5.12. 
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5.5,19 Subroutine DERIV 

This subroutine computes the total derivatives and the partial 
derivatives g. ^ described in the analysis presented in Section 3 of Reference 
3 * 

The implicit integration method used to integrate the differential 
equations governing the chemical system, i.e. 

y i = V x,y l* ' ' * y NSP+3* i=1 ' ' * " NSP+3 

v/here the variables . 


y t i=l, . . NSP+3 

are V, p, T, i=l, . . NSP 

respectively, requires evaluation of the Jacobian of the system, i. e. 


Bf 


p ij” ay 


i 


i=l, . . ., NSP+3 
j=l, . . NSP+3 


Subroutine DERIV computes only certain of the p (those taken 
with respect to C^) and the others are computed in subroutine FLU. 

Also calculated by DERIV are the reaction rates, kj, and the net 
production rates, X^. 

The generalized chemical reaction which is handled by this sub- 
routine is defined by: 

NSP NSP 

. ,?! nj V 

where Mj represents the i th chemical species. 

The reverse reaction rate constant is defined by the equation: 


kj SK(J) 


=a J • T ~ n .j ■ exp (-bj/T) 


The net production rate for a reaction is given by: 

C i 


T NSP NSP V 'l 

X. Xffl -[ K l- W '->■ -H, C i “J 
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where: X depends on the order of the reaction* 



NSP 


with 

M J = 1=1 mM i.i ' C 1 

for reactions requiring a third body 

and 

M j = 1 

for all other reactions 


The net individual species production rate is given by the equation: 
dC. _ ^ 

-sr wo -v fa 


where: 

K a = (Mw 1 • p • r*)/V 

*ij E V v u 

The partial derivatives of the net species production rate with 
respect to: the chemical species; the gas velocity; the gas density; and the * 
gas temperature are: . 

_ T »x, 

BT(I,K) = Kj • i = l NSP 

1 k = 1, . . NSP 

1 dC i 

PHI(1, 1) = -“7"-^- i = l NSP 

i dC. _ y 3X 

PHKI.2) = ■ + j=i ^ 1 = 1 NSP 

l 

_ y ax 

PHI(1, 3) = K i j=l ST i NSP 


* X = (Vjj - v^) so that X = 0 for binary exchange, X = 1 for most 


i=l 

dissociation recombination reactions. 
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The subscript notation used above is: 

1 *= Species subscript 

j « Reaction subscript 

i = Total number of chemical reactions 

m * Number of reactions requiring third body terms 

NSP = Total number of gaseous species 

Subroutine DERIV calls subroutine FLU to* calculate the derivatives and 
partial derivatives of V, P, and T. In the event of solidification , this routine 
also recalculates the condensed phase species derivatives and their partial de- 
rivatives and recalls subroutine FLU to recalculate the flow derivatives. See 
Volume I, Section for the description of this procedure and the equations 

calculated. 
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5,5.20 Subroutine FLU 


This subroutine computes the total derivatives and the partial 
derivatives a 4 and 3^ for the fluid dynamic equations. While the flow is 
subsonic, pressure defined fluid dynamic equations are used. When the 
flow becomes supersonic, area defined fluid dynamic equations are used. 

The summation terms, energy exchange term B, the diabatic heat addition term 
A, the Mach number, and all the partial derivatives of these terms are com- 
puted. For a subsonic integration the pressure and its derivatives are obtained 
from the subsonic pressure table. For a supersonic integration, the area ratio 
and its derivatives are computed from the input geometric constraints. 

The calculations logically fall into three types: a) Those done for 
all integrations; b) Those done only for subsonic integration; c) Those done 
only for supersonic integration. The following will adhere as closely as 
possible to a sequential description of the computations. 

The operators $ (i , 1) , 1 = 1,2,3 are defined as 
$ (i 1 1) « P (C., V) 

ft(i,2) = )3 (C A , p) 

$(i.3) = fi (Cj.T) 

dC { 

The total derivatives, f i = -gj- , for i = 1, . . . , n are computed as 
tt> 4 r* 


where n is the total number of species, NSP. 


Computation of the Summation Terms and their derivatives; 
First Summation 


81 

SI 

- 

1 11 C 
I -S - 
R 1=1 

dSl 

dV 

DS1V 

a 

1 11 

B. * 

dSl 

dp 

DS1R0 

St 

rli *i 

dSl 

dT 

DS1T 

a 

1 11 

i Si • , 

dSl 

dCi 

- DS1C{I) 

- 

i . fs 
r Lj=i 

Second Summation 






, n 

82 

02 

B 

-L-s 
r-t i=: 

dS2 

dV 

DS2V 

S 

1 11 

R-T i=I 

dS2 

dp 

DS2R0 

= 

i n 

r • t i=: 

dS2 

dT 

DS2T 

a 

R*t£1 


'1 


R. 


ffr DS2C(I) 


R •[|5l , ’(C ) .0 1 ) V S1 ' K l] " 


dc. 


1 £ *- , — 1^ 1 S2 

R-TiM t_*(1.3) ' h l + dx Cp l j T 


1 


^ n. ^(C ,C )' h l 

R-Lf-l I » 


*R is the gas constant/molocular wt. of species i except for condensed phases 
whfcre R^O. 
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Computation of the Energy Exchange Term B and its Derivatives: 


B 

BB 

- 

1=1 . 

7 

S2 


dB 

dv 

DBBY 

- 

1=1 . 

7 

as2 

av 


dB 

dp 

DBBR0 

MS 

1=1 . 

7 

as2, 

ap 


dB 

dT 

DBBT 

= 

1=1 . 

7 

as2 

ai 

. J2 . & 

y 2 a T 

dB 

dCi 

DBBC(I) 

= 

1=1 . 
7 

as2 

ac t 

+ S2 . a* 
y 2 dc i 

Computation of the Diabatic 

Heat Addition Term A 

and its Derivatives 

A 

AA 

SET 

Sl-B 



dA 

DAAV 


asi 

dB 


dV 

SS 

av 

dV 

• 

dA 

dp 

DAAR0 

= 

asi _ 
8 p 

dB. 

dp 


dA 

dT 

DAAT 


asi 

dB 


SS 

St 

dT 

‘ 

dA 

d^ 

DAAC(I) 

= 

asi 

3 c i " 

dB 

dc i 

i= 


I 
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' 3 








m 


mammae* 


Computation of the Mach number and its derivatives: 


M 2 

XM2 

II 

»-3 




*M 2 

av 

DM2V 

_ 2 -M 2 

V 




3m 2 

3T 

DM2T 

2 2 

Ml _ MT 

T y 

±7 

3T 



3M 2 

* C | 

DM2C© 

XK i r ay i 
- - M Lac; r 

R in 

+ TJ 

i ** 1, . , 

> • , n 


For the subsonic portion of the nozzle, pressure defined fluid dynamic 
equations are used. The pressure, and its first and second derivatives are 
computed via interpolation in the pressure table generated by subroutine PRES. 
The Subsonic Gas Velocity derivatives are computed: 


dV 

dx 

FNX(l) 

1 

sr — ■ • 

P‘V 

dP 

dx 

d[FNX(l)] 

ax 

AL(1) 

1 

“ " 'p V * 

d 2 P 

dx 2 

0(V.V) 

BETA (1,1) 

1_ 

V 

dV 

dx 

U(V,p) 

BETA(1,2) 

i 

p 

dV 

dx 


The Subsonic Gas Density derivatives are computed: 




dx 

afrNXfe) ] 

dx 


FNX(2) 

AL(2) 
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fi{ P.V) BETA(2 , 1) p • av 

/J{p,p) BETA(2,2) ■- ^ ~ p ' bf 

P(p,T) BETA (2, 3) =- P • f* ~ p^2‘ fj * £ 

% 

Wp.Cj) BETA (2 , 1+3) = - “ 2 ^ ' ^ * JJ- p-|*' 1=1.. ...n 
The Subsonic Gas Temperature derivatives are computed: 


dT 

dx 

FNX( 3 ) 

T • 

SlMffl 

dx 

AL( 3 ) 

ZzL . 
y 

/KT.V) 

BETA( 3 , 1 ) = - 

T • 

0 (T.p) 

BETA( 3 , 2 ) = - 

T • 

/J(T.T) 

BETA( 3 , 3 ) = 

I . dT 
T ’ dx 

^(T.Cp 

BETA( 3 , 1 + 3 )= 

T • . 


~Z± . i . dP _ -> 
L y P dx J 


1 !" d 2 P /flP N 2 I 1 

p ’ Ldx 2 

dj 

dv 


dB 

dp 


+ T 


_i_ d£ ? 2 L _ T 

y 2. p dx 8T 


dB 

dT 


r JL dp M 

L y 2 . p dx dC A - dc^ 
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For the supersonic portion of the nozzle, area defined fluid dynamic 
equations are used. The area ratio, and Its derivatives are computed according 
to the Input geometric constraints. 

Area ratio and Its derivatives: 

1) On the circular arc of radius (input Item RWTD) defining 
the downstream throat region, X s ^ tangen t 



2) For a conical nozzle and X> X 


tangent 


da 

dx 


j't + ( x " x c) tan °t J 

. 2 + ^x-x^ tan °t j' tan ®t 


4 " 2 tan^ 0 

dx* S 
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3) For contoured, parabolic, and circular arc nozzle and x> ^ an g en t 


- Y 2 



where Y, dY/dx are computed via interpolation in the tables of Y and its deriva- 
tives generated by exact calculations or by subroutine SLP. 

The Supersonic Gas Velocity derivatives are computed: 


dV 

dx 

FNX(l) 

JL . ['I da _ A “1 
m 2 -i La J 


ilm xjiJl 

dx 

AL(1) 

V_ . i . rA _l.p da. 
M 2 -l a L dx 2 8 V * 

■)’] 

/»(V,V) 

BETA(l.l) = 

I . dV JL . dV . dM? 

V dx m 2_ : ‘ dx dV 




V . dA 



■ 

2 by 

M -1 


W. p) 

BETA(1 , 2) = 

V dA 

"m 2 -i’ ap 



BETA(1 ,3) = 

1 dV . dM 2 V . 

" m 2 -i ■ ** aT m 2 -i 

dA_ 

dT 


BETA(1 ,1+3) = 

1 . dV . dM 2 JL . 

“ m 2 -i ac i ’ m 2 -i 

dA 

dCj W " 
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Ihe Supersonic Gas Density derivatives are computed: 

. FKX(2) » - p{ ” A ) + A ] 


MFjTO J ^(2) 


M 

p • -y~ 

M -1 


1 fd 2 a 1 ( ^da "fl 
a ’ L dx 2 ' 8 


Pi p,v) 

BETA(2,1) = 

r 1 f 1 da A 3M. _J_ . \h. ~j 

> V.,l> ■ = -V *> M z -1 wJ 

Pi p.p) 

BETA(2 / 2) = 

JL.de. + _fi_ . |A 
p M z -1 p 

Pi P.D ■ 

BETA(2,3) = 

r 1 ,M da y&M 2 1 • — 1 

P<L (m 2 -i ) 2 dx J dT m 2 -i 9tJ 

fiip.GJ 

BSTA(2,l+3) = 

- 1 /■ 1 da -\dM 2 , _i .^A ~j 

11 ' i- (M 2 -!) 2 ^ * dX 790 ‘ mW 



l r l / • • • / 


The Supersonic Gas Temperature derivatives are computed: 

FNX{3) - - T • [(r-D • M T' K a fc ” k J + B j 
dx L M -1 


mmi alO, 


S(T,V) BETA(3 , 1) = T 


0(T,p) BETA(3 , 2) 


2 , _ ,2 

M y-1 ,! d_a. 

M 2 -l 8 dx 2 

_I 'da Y“| 
a V dx J J 



r y - 1 V 1 da A 

»v * r V-l 

. Ik 

a_B*i 

V-D 2Vadx ' 

dV 

“av J 

, r i m 2 .M 

aj i 



L y ‘ 1 2, dp 

L M -1 

ap j 
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P(T,T) BETA(3,3) 


I.M 7 . 1*4 . — lli. 




r y-i r Ida -s sm 2 
L^zUdx " J 3T 


+r 


, , M Z 3A 93 
-1 • — r* • »r - sr 


M 2 rl 


9T 91 


/3(T, Cj) BETA P , 1+3)= 


_ mL . ri da 

„.2 \a dx 
(M -1) 



T.r * 

L (M 2 -!) 2 ^ “ * 



M 2 , 3A 


9B f\ da a >iZ_1 


1 = 1, 


• • • / D 
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v: 


5 >5 .21 Subroutine GTF 

This subroutine computes the effective gas constant, gaseous heat 
capacity, y, dy/dT, dy/& from the following formulae: 

NSP 

R - S 0,-R, 

1”1 

NSP 

Cp - E C.Cp. 

i»l 1 1 . 


221 


£J2_ 

Cp-R 

2L.fr.-U- 

CP 


NSP 
E C 
1«1 


acpj 

l dT 


221 

*C 


y. fr-1) 


1 


. r~ 

Lr 


Cp 


11 
Cp j 


1=1, . . n 


For condensed phases, R^O , to account for the assumption that the particles 
exert no pressure on the gas. 
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5.5.22 Subroutine IAUX (1IL, H, QK, RK, IX) 

This subroutine performs implicit integration according to the method 
discussed in Section of Reference 3. The Increments for the chemical species 
concentrations and the fluid dynamic variables at thelorward point are cal- 
culated by solving the appropriate implicit finite difference formulas. Subroutine 
IAUX also performs explicit integration, using a modified Euler method, when 
the gas temperature falls below an input value. 

The calling sequence parameters are: 

HL - last integration step size 

H - current integration step size 

QK - last increments for variables 

RK - computed increments for variables 

JK - 1 initial 3 steps 

2 general step 

3 special step 

4 restart step 

The total derivatives, f. , and partial derivatives, p. . at the 

i,n i , J , n 

back point are calculated in subroutines DERIV and FLU. 

The special step calculation is used at print stations, in halving 
the step size if required, or for integrating to specific calculation stations. If 
the special step calculation is used to determine the properties at a print station, 
the calculation is resumed using the general step calculation and the previous 
step size. 

After each integration step, subroutine IAUX obtains the derivatives 
at the then current axial position. 
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For Implicit Integration the equations used are: 
Initial Step and Restart 




°i,o h 


N 

+ r 


e. . K 

i.j.o j 




General Step 

^i f n+l 

Special Step 


ik.* 2 1 


f . + a. h 

i,n i,n 


N 

+ r 


i,j,n J ,n+l 



n+1 


i,n+i (2h +h )-h 
ui'JL n n 


N 


k . + 

i,n 


f . + a. h 

i,n i t n n+1 


+ J-1 ’ h n+ l ( 


h..+h ) 
n+1 n 


For explicit integration the above equations are used deleting the 
partial derivative terms a and £ . 

If the option to generate input tables for the Turbulent Boundary Layer 

Nozzle Analysis Computer Program was selected, tables of M, P/P T/T , C , 

c c p 

V, p, are tabulated using subroutine TABGEN. While the code remains in this 
routine to do the above, the SPP code does not use these values in any way. 
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5.5,23 Subroutine INT 


Provides control for the implicit Integration procedure , determines 
the proper set of nonhomogeneous equations to solve, and, after each 
integration step, computes the next integration step size according to 
the following relations: 


w - 2 W 


h n+2 * 2 n+1* 



k t,n41 ~ 2k j.n 4 k l.n-l 
3k, .... - k. 


i.nfl i,n 


< 

MAX 


A»n+1 i.n l.n-1 


3k - k 
1 , 1*1 K i,n 


MAX 




k - 2k, + k. , 

_ i,n i.n-1 


3k - k, 
i,n+l i,n 


MAX 


6 _ 

10 


6 


6 


On option, (JF=1) only the fluid dynamic variables are used in 
determining the next integration step size. 

If the step size is halved for the fourth step, the integration is re- 
started using one-half the original step size. 

The correspondence between equation number and physical property 
is: 

Equation Number Property 

1 Velocity of Gas 

2 Density of Gas 

3 Temperature of Gas 

4 *♦ NSP+3 Gaseous species mass fraction 

(1 -*NSP) corresponds to (4 -♦ NSP + 3) 


When the flow is supersonic, continuity is used to control the integration step 
size to insure that: 


C pVA) " ( pVA)^ 
<-PVA ) n+1 


< CtfNDEL 


where CGfrJDEL Is an Input relative criterion. 
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5 , 5 >24 Subroutine LESK (Y) 

This subroutine Is a single precision linear equation solver which is 
used to perform the matrix inversions required by subroutine IAUX. Gaussian 
elimination is used with row interchange taking place to position maximum 
pivot elements after the rows are initially scaled. 


l 

s 


l 

T 
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5,5.25 Subroutine MAINin 


This subroutine provides the overall logic control for the one-dimen- 
sional kinetic expansion. The following functions are controlled: 

1) Variable initialization 

2 ) Option to start the kinetic expansion from equilibrium 
throat conditions 

3) Controls of the integration to hit specific* area ratios, the nozzle 
throat point, the nozzle tangent point, and the requested end point 

4 ) Controls of the switch from the subsonic pressure^iefined equations 
to the supersonic area defined equations when (M ^1.02) 

5) Controls the switch from implicit to explicit integration. 

For the normal mode of operation of the program, this subroutine locates 
the throat in the following manner: 

The gaseous mass flow per unit area (pv) is calculated and stored as a 


: 


* 

* 


i 

4 


function of nozzle axial location for the present and past integration step. When 

(pv) n + l < (pv) n 

where n refers to the n* h integration step, the throat location is calculated from: 


(X - x ) 2 - 
n n-i 

(p V ) n+ i -(pv)J 

+ ( VrV 2 - 

(pv) -(pv) 

n n-i. 

2 • 

’ (x n + rV 

• r (pv) -(pv) 

l n 1 

l-l 

"(X ~X 1> 
n n-1 

•[(P v kr<r.v) n j 


and the n+l** 1 Integration step is repeated using a step size of X* - X r to deter- 
mine the throat conditions . 

T6 prevent the location of a false throat due to roughness of an input 
pressure table, ten integration steps are required before the throat will be 
sought. 


Through the downstream throat radius of curvature the step size is 
controlled so as to be less than or equal to RWTD* SIN (THETA) /25.0 f 
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Subroutine MAIN1D also contains the logic to control the integration through 
solidification of multiple condensed liquid phases if they are present. This logic 
allows the integration procedure to hit the beginning and ending of solidification 
exactly and to turn on the solidification equations in subroutine DERIV via the flag 
IFMELT in COMMON LKMELT. 
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5. 5 >26 Subroutine OUTPUT 

This subroutine provides conversion from internal computational 
units to output engineering units and the calculation of performance para- 
meters. The following output parameters are computed by this subroutine: 

*« 

The pressure (in PSIA) is computed from: 

P (PSIA) = p /«33.056 

The gaseous species mole fractions are computed from: 

R 1 

c l.» ■ T • c l 

The gas molecular weight is computed from: 

Mw = 49721. 011/lR 

The percentage mass fraction change is computed from: 

n 

% A (Mass Fraction) = 100.0 • (1.0-2 C.) 

i=l 1 

The gas heat capacity is computed from: 

Cp (BTU/LB-°R) =3 . 9969 • 10 -5 • Cp 
g y 

The gas enthalpy is computed from: 

n 

H (BTU/LB) =3 .9969 • 10 -5 ' 2 C, • h. 

8 1=1 1 1 

At the throat, the characteristic exhaust velocity (ft/sec) is com- 
puted from: 

c * = A_ 

w p* • V* 
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Xhe vacuum specific impulse Is computed from: 


ISP. 


V + p/p - c * 

c _ g- 32.174 


VAC 9 

Ihe vacuum thrust coefficient Is computed from: 


I 


r VAC 


s Pvac 

~~C*~ 


Ihe percentage enthalpy change Is computed from. 


%ah t = 


100' (HREF -HREF) 

c 


V 2 /2 


where 


NSP 


HREF * p l C A - ^ + V /2 


HREF is HREF evaluated at the initial condition for the 0DK 
Integration (i.e. at the initial contraction ratio, ECRAT). 

In addition the actual molecular weight of the mixture of gases and con- 
densed phases is also printed out. 
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5,5.27 Subroutine PRNTCK 

For the option to print starting at step ND1 # printing every ND3 rc * 
step up to step ND2, this subroutine checks whether or not the current step 
should be printed. If it is to be printed this subroutine calls subroutine 
0UTPUT. 
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5.5.28 Subroutine TAB GEN (IFLAG, L TABLE, XTAB, YTAB. LUSED, X, Y, IERR0R, N 

This subroutine records a tabular function (X,Y(NY)) in tables of fixed 
length. The first and last event will always be tabulated and the table will 
either contain all of the values specified or will be at least half full. Once 
the number of events exceed the table length, the table will be repacked by 
the deletion of every other table entry and tabulation will proceed choosing 
every 2— ^ event (N=0, 1,2... etc. , where N is the number of times the 
table must be repacked). The table spacing will be a power of two except for 
the last event which will always be tabulated. 

The calling sequence parameters are: 

IFLAG - denotes type of entry to subroutine 
■ - 1 first entry 
* 0 normal entry 

*= + 1 last entry 

LTABLE - length of tables available for tabulation 
XTAB - table for tabulation of the variable X 

YTAB - table for tabulation of the variable 

LUSED - number of table entries currently used (output) 

X the variable X 

Y - the variable Y 

IERR0R - error flag 

NY - number of Y variables to be tabulated 

Note: One Dimensional Mach Number Tabulation Procedure 

At the initial axial position, X and Mach number are recorded. 

TABGEN is then used with LTABLE=50 (assuring 25 saved values). The last 
recorded values are the end values for the transonic tables. 
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5,6 Link 50 Subroutines 

This overlay and subroutines contained within are the essence of the TD2P 
module. The majority of this section is from Reference 4 with the major modifi- 
cations restricted to Sections 5,6,2 through 5,6,12. However, almost all of 
the subroutines have been modified to some extent in order to bring this calculational 
module upto date. That is, this calculational module has been updated from 
F0RTRAN II to F0RTRAN IV by the extensive use of labeled common blocks and the 
reductions in the amounts of equivalences to blank common areas. 

Other modifications made to the TD2P module include the extensive trans- 
mittal of data to reduce the amount and complexity of data needed to execute this 
module. 


5,6,1 Program LINK5Q 

This subprogram is mainly intended to facilitate the conversion of the 
SPP code overlay structure to other computer systems. The only function of this 
routine is to call subroutine TD2P, 
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5,6.2 Subroutine DIST 


Input: 

Output: 


D, N, <j g 

u \ * L 

J m T 


M, ... N 


Description : 

Given an arbitrarily large number of spherical drops of mass mean 
diameter, D, which are assumed to possess a logarithmic normal distribution 
about D, then the mass density is given by: 



d D 


[(2 irl /2 D In <r g ] 


exp 




2 


where 

D is particle diameter 

m,p is the total particle mass 

< 7 g is the geometric standard deviation 


It follows that 



The purpose of this subroutine is to determine a set of N drop diameters 


j=l, ... N 

such that 



The method described above is the same as presented in References 10 and 11. 
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5.6.3 Function EISP 


Function calculates the error function 

EISP = c - c (x) 

used in finding the one ^mensional - zero lag I as a function of area ratio, e . 

sp 

The independent variable, x, is either the gas temperature (before or after solidi- 
fication of the condensed phase) or the fractional amount of solidified condensed 

phase (during solidification). The I corresponding to the independent variable 

sp 

is also calculated. 

This function is called by subroutine SEEK which is called by subroutine 

ISPID. 
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5,6.4 Function FAM 


This function calculates the error function 
FAM (M) = f - C ( M) 
where c = requested area ratio 

€ (M)= area ratio at the Mach number, M, given as the argument of function 
FAM. 

The FAM function is called by subroutine SEEK (Section 5.6.8) which is called 
by TBLEDG (Section 5.6.9) to determine the subsonic-transonic Mach numbers 
as a function of area ratio for communication to the TBL module. 

The error function is calculated using the following: 

A*/A = ASA = *±i M (1 + x y ~ Ms) 2 ^V _1 )/(v +1 ) 

FAM = e - 1.0/ASA 
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5,6,5 Subroutine FIND 


The purpose of this subroutine is to perform linear or quadratic 
interpolation from a table. 

Caning Sequence 
CALL FIND (N,T,X,Y) 
where 

N is a flag 6uch that: 

N • 0 implies linear interpolation 
N jt 0 implies quadratic interpolation 

T is a table of the format specified below 

X is the argument 

Y is the result 

Usage 

The table, T, must be a column of a FORTRAN array. The format of the table 
is: 


T(l) a temporary cell used by FIND 

must be 0. initially. 

T(2) minimum argument 

T( 3) function of minimum argument 


T(2n) * nth (maximum) argument 

T(2n+l) function of the nth argument 

T(2n+2) 0. signifying the end of table 

T(2n+3) 0. signifying the end of table 


If FIND is entered with an argument above the table maximum or below the table 
minimum, it will store the function of the maximum or minimum argument, respect- 
ively. If FIND is entered with an argument which is between the first or last 
two table arguments, linear interpolation will be used. 
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Functional Description 


This subroutine saves its place In a table and > on an option, either inter- 
polates by the linear interpolation formula: 

y . (*-*!> +*! * 1 <X **2 

*1 

or performs a four point quadratic interpolation by the formula: 

2 

y - Afx-Xj) ♦ ♦ C x Q <x x < x s Xg < x^ 

where 

d - x 2 -x 1 

* - V*0 
f mx r*i 

a - (y 2 -yj) ( e -U/* * + y 3 

d(e-f) + e 2 ♦ ? 

B - (yjj-y^/d - Ad 

C - y x 

The above quadratic formula Is constrained such that the function y{x) passes 
through points (x^, y^) and (x,, y^) and has the property 

[y(x o )-y 0 ] - [y(x 3 ) - y 3 ] ■ o 


Note: There is, depending on the computer system used, a possible entry point 
name conflict between this routine and the subroutine of Section 5.1.6 by the 
same name. These two routines are not the same. See Appendix A for conversion 

aides. 
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5,6.6 Subroutine ISPI D 


This subroutine calculates ; for a given area ratio assuming one dimen- 
sional equivalent perfect gas flow with zero velocity and temperature lags between 

the gas and condensed phases. This I is referred to as I in most of the 

S P sp lDO 

documentation concerning the SPP code. # 


The calling sequence to ISPID is CALL ISPID (E, ISP , XI) 

where E = area ratio, c , at which The^I is desired 

sp iDO 

ISP = I on exit from ISPID f 

Sp lDO, 

XI = initial guess at the independent variable 


Before and after solidification of the condensed phase, the gas temperature 
is used as the independent variable to iterate on the solution, however, during 
solidification, the fractional amount of solid to total condensed phase is used. 

This routine utilizes subroutines SEEK (Section 5.6.8) and EISP (Section 
5.6.3) to determine the root of the appropriate equations. 

The constants used in these routines are calculated in subroutine AGP 
(Section 5.6.12). 


B 
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5.6.7 Subroutine SEEK 


This subroutine utilizes subroutine ITER to find the root of the equation: 

f 00 = o 


Subroutine ITER predicts new guesses of the root, x^, given old values of the error, 
f Q , and guess at the root, x Q , using the method of secants. The SEEK subroutine 
' attempts to take maximum advantage of the properties of the secant method. (See 
the write up on subroutine ITER). 


Calling Sequence 


where 

XANS 

X0 

MAXIT 

EPSI 

FCALC 

IB 


B(I) 


Call SEEK (XANS, X0, MAXIT, EPSI, FCALC, IB, B, IT, IERR) 


the root of f(x)=0 or the last guess for x if the 
maximum number of iterations has been exceeded. 

initial guess at the value of the root. 

maximum number of iterations allowed in finding 
the root. 

a so^tion to f(x)=0 is claimed if 
Jf (x)j < EPSI or 

I Ax | < EPSI/100 . 

an external function subprogram which calculates 
the error, f(x). The calling sequence is F=FCALC(x). 

a bounding flag such that if 

IB = 0 no bounding. 

IB = 1 the root is assumed to lie within the 
values of B (1 ) and B(2). 

IB = -1 the maximum absolute change allowed 
between successive guesses of the root 
is constrained to be no larger than 
|B(l)*|x o | + B(2)| 

a vector of length 2. The use of this vector is de- 
scribed above. 


IT = . the number of iterations required to obtain a solution. 

IERR = a flag such that if 

IERR=0 a solution was obtained. 


IERR=I failed to converge on the root in MAXIT 
iterations. 
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5.G.8 Subroutine TBIEDG 


This subroutine calculates the subsonic-transonic edge conditions for the 
trubulent boundary layer module, TBL. The number of points, NTBL, and the initial 
location of the subsonic start conditions, XITBL, are determined by user input in 
the $TD2 namelist (Volume in. Section 2.9). 


The method used in constructing the edge tables for TBL is as follows: 

1 . The region from XITBL to the MOC Initial line is broken up into NTBL 
equally spaced points. 

2. Subroutine SEEK and function FAM are used to calculate the 1-D equi- 
valent perfect gas Mach number, M _ at the area ratio where the 
initial line intersects the wall, 

K U 

3. The axial position, Z , in the nozzle which corresponds to an inlet 
area ratio of nine is ftcated. 


4. The edge Mach numbers are then calculated from 


M = M in 
e ID 


il + 


Z-Z 

o 

zpr 

i i o 


M. -M ini 

■ A t IPl j 
M 1DU 


I 


where the subscripts 


i i = initial line 

ID = 1 dimensional Mach numbers calculated using SEEK and FAM 
as a function of axial position, Z. 


In addition to the above, subroutine TBLEDG also recalculates the stagna- 
tions pressure and temperature communicated to the TBL module. The values of 
P Q and T q are calculated from the static P and T at the initial line wall point in 
order that edge conditions calculated by TBL would more closely resemble the 
values calculated by TD2P. 
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5.6.9 Subroutine TD2P 

Subroutine TD2P is the master subroutine of the TD2P module. It controls 
the initialization of variables, stores in the proper locations values transmitted 
from other modules, reads in user input in the $TD2 namelist, and directs the flow 
of the computation. 

In the case that particle sizes are not input to the module, TD2P calculates 
the mean particle size (in microns) as 

d = XK • PC PEXP • (») XIEXP (i _ EXP (XL *XL STAR)) 

P 

*(1 + 24 • XD • r*) 

£ = 100 (WPWGT/((1 +WPWGT) • XMLW)) 

where XK, PEXP, XIEXP, XL, and XD are constants set in DATA statements or 
they may be input in the $TD2 namelist. 

PC = chamber pressure (psia) 

XL STAR = L* (in) 

£ = moles of condensed phase/lOOgm. 

WPWGT = m gas/m condensed phase 

XMLW - molecular weight of the condensed phase 

The particle size distribution is then calculated by calling subroutine 

DIST. 

Overlay 5, 1 (LINK51) is called to calculate the average gas properties 
used by the rest of the TD2P routines (see subroutine AGP, Section 5.6.12). 

The following items are then computed and made available: 

P 



5-113 





Overlay 5,2 (LINK52) Is then called for execution of the Axlsymmetric Sub- 
sonic-Transonic Subprogram and then the following items are computed and made 
available: 


C 

c 




R) 


N 




P 

r 


m 

P 



(T R) N 
g o 


r* 


r 2 

p j 


j-=i, — , j 


max 


after which Overlay 5,3 (LINK53) is called for execution of the Axlsymmetric 
Supersonic Method of Characteristics Subprogram. 


- 




i 


i 
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5,6,10 Program LINKS 1 


This subprogram Is mainly intended to facilitate the conversion of the SPP 
code overlay structure to the computer systems. The only function of this routine 
is to call subroutine AGP, 
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mm 




5 . 6 . 1 1 Subroutine AGP 


The purpose of this subroutine is to calculate average gas properties based 
on the results of a one-dimensional calculation using variable thermodynamic pro- 
perties (usually , the results of an 0DE calculation). These properties are then 

used in the calculation of one-dimensional and two dimensional I in the rest 

sp 

of the TD2P module. The method used is outlined below. 


Using Newton's method, an expansion coefficient, v , is found such that 

a polytropic expansion from the chamber temperature, T gQ , to the solidification 

temperature, T , will occur at the input expansion ratio, e . The iteration equa- 
pm *pm 

tions used are 


where 


e 


-(1+1) _ -(i) / c c pm \ /dr 
v - V ( - € } €/ dy 


'<f - » 

v 


1/2 




-X.+-L 

2(y-D 


and 


It 

d y 


(V- D 2 




Next the value of the specific heat, C , for the gas phase is calculated 
from the one dimensional energy equation as 




2 V V 



i.e. the reference enthalpy at stagnation is taken as zero. The flow velocity 

at the T is given. It is also assumed that C . for T ^ T _ is also known, 
pm 3 pi pm 

An expansion coefficient for the gas phase is next calculated assuming 


no lag: 


1 - (y - 1) c Pi 


iat C 
9 pg 


so that 


Next the value of C^, the particle specific heat after solidifi- 
cation (T < Tp m ) # is calcuiated from the energy-equation as 


w (1* - T ) 

p' pm A ps' 


i 1 + 


u - c fr - t ) 

gs w pg ' go ps' 


C pl “ V - -- 

w g 


The values of T and U are input values corresponding to some point 
ps gs 

after solidification. The above value of C ps is then checked against the actual 

value of C at the solidification temperature. If the calculated value is more 
ps 

than 10% different than the actual value, the actual value is used. The reasons 
for the above test are 

a) it prevents ridiculous values of C from being used 

ps 

b) it prevents the results of the TD2? calculation from being too sensi- 

tive to the area ratio selected at which values of T „ and U „ are ob- 
tained ps gs 

If the Prandtl number, Pr, is not input or transmitted from the 0DE module, 
it is calculated, using kinetic theory, as 

Pr = 4 V /(9 V - 5) 
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5.6.12 Program LINK52 


This subprogram is mainly intended to facilitate the conversion of the SPP 
code overlay structure to other computer systems. The only function of this rou- 
tine is to call subroutine PARTIL. 


5-119 


5,6.13 Subroutine ABCALC 


Coefficients are calculated for first and second order transonic 
approximations : 


First Order; 
C - 


u 2 - 1 
0 


- k-1 2 

1 ~ TTT U 
k+1 O 

C a 
o 


20 AC. 


u a 
o 


21 2C, 


V21 0 


40 8C, 


Second Order; 


A 201 

=5 

, k-l 

2 i5r u o a 




4A 201 a 20 + C o S 


b 21 

VC 

4C 1 


A 102 

B 

2 

a + u 3 

0 


* 


2u o a6 + A 102 a + 

4A 201 a 21 

b 22 

m 

4C 1 


A 120 

b 

2u b 01 

0 21 


A 220 


0 k-l 

2 k+T u o a 21 


A 320 


k+r u 0 a 2 o 




2C o b 22 + A 120 a 

+ 4A 220 a 20 + 2A 320 a 21 

b 40 



16C X 

A 121 

=3 

2a 21 “ + 4u o b 22 






*321 “ k+1 u o a 21 

2u o a 21 8 + 4u o b 22° + A 121 a + 4A 220 a 21 + 16A 201 a 40 + 2A 321 a 21 
b 41 “ 16C 1 

A 140 “ a 21 2+2u o b 41 
A 340 “ k+f U o a 40 





jA 


5.6414 Subroutine CCALC 


Coefficients are calculated for the third order transonic approxi- 
mations: 


Third Order; 


B 

k-1 A 


°202 

k+1 A 102 


c 

C o y/2 + 4A 201 b 21 

+ 4B 202 a 20 

C 22 

4C 1 


B 103 

■ u i + aB 

0 j 


c 

4A 201 b 22 + 4B 202 3 

21 + A 102 6 + B 103 a + u o ay 

C 23 


4C, 


n , k-1 2 

B 120 " 4 k+1 a 20 


t. k-1 

*220 " k+1 A 120 


*40 


B 12Q a + 4B 220 a 20 * 2A 320 b 21 * 2C o C 22 
16C, 


B 121 ” 4u o C 22 + 2b 2l“ + 8 k+1 a 20 a 21 


r k-1 , 

°221 k+1 121 


B 321 ■ iST (a 20“ + u o b 21> 


'41 


16C. (4A 320 b 22 + 2A 321 b 21 + 16A 201 b 40 + 4A 220 b 21 


+ B 121 a + 4B 221 a 2Q + 4B 22Q + 2B 321 a 21 + A^B 


+ 6C o C 23 + 4u o C 22 c) 


"122 


6u o C 23 + a 21 6 + 4b 22° + 4 M a 2l' 






(a 01 a + u b 00 ) 


*42 


B 


"80 


16C X (B 122° * 4B 221 a 21 " 16B 202 a 40 * 2B 322 a 2* * A 121 6 


+ 2A 102 b 22 + 16A 201 b 41 * 4A 220 b 22 * 4A 321 b 2? 


12u o c 23 a + « o a n y) 


k-l 


140 


2u o C 41 + 2a 21 b 21 + 16 k+I a 20 a 40 


k W . 

a 240 * k+1 A 140 


B 340 * k+1 (s 20 a 21 + 2u o b 40 } 


C 60 " 36C 1 (B 140 a + 4B 240 a 20 15B 220 a 40 * 2B 340 a 21 + 2A 12Q b 22 


+ 16A 220 b A o + 2A 340 b 21 * 4A 320 b 4i, + 2C o C 42 + 4u o a 21 C 22* 


k-l 


141 


4a 21 b 22 + 2b 41 a + 16 k+1 a 21 a 40 * 4u o C 42 


B 34l " k+1 (a 21 + 2a 40° + 2u o D 41 ) 


JL , 


+ ^ A ^'>i^ai + ^6A* n b- n * * 


“61 36C X v "34C U 22 321^41 ^ w 201 60 >20 43 


+ fc 14l a + 4B 240 a 2i + l6E 221 a 40 + 2B 341 a 21 * A 140 £ 


2A 121 b 22 + 4u o C 42 a + l2u z a 21 C 23> 


160 


2a C,.. + 2a 01 b. , + 16 T7~ a, n 2 
o 61 21 41 k+1 


B 360 • k?T (2a 21 a 40 + 3u o b 60 } 


64C 1 v 16.T 


+ A6B 240 a 4C + 23 360 a 21 + 2A 140 b 22 + 36A >™ b ‘ 


220 60 


~ 4A 340 D 41 * 4u o a 21 C 42 ) 
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5,6.15 Subroutine DCALC 


A dummy subroutine containing only a RETURN statement. This sub- 
routine is reserved for the calculation of coefficients for the fourth 
order transonic approximations. 
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5.6.16 Subroutine FCALC 


This subroutine calculates the boundary condition functions, f 
from coefficients and derivatives determined by subroutines ABCALC, 
CCALC, and PCALC: 


u - u 


t dr 

f 0 * u — - v 
2 o dx 

2 

f du dr , d r _ dv 

*3 dx dx U o , 2 ” dx 

dx 

d 2 u dr - du d 2 r 

t . * — - — + 2 r + u « 

4 . 2 dx dx , 2 o,3 

dx dx dx 


.3 ,2 

dr d v 


dx 4 


d^u dr . d 2 u d 2 r _ du d^r d^r d\ 

, 3 dx + 3 . 2 , 2 3 dx , 3 u o , 4 " . 3 

dx . dx ax dx dx dx 


If a special throat expansion (see Step 3, subroutine PARTIL) is 


being performed, the following variables are computed: 


* 0 o , o 

u ■ u +a x +8 — + yt~ 
o o 2 '6 

2 

x 


<** - a + fix +Y7 2 
o L 



If transonic conditions in the nozzle inlet are to be fared to the 

throat conditions, the following variables are computed: 

2 3 

. , . (z +x ) . (z +x ) 

A k fJk y n k 'wo 

u *= u +• a (z +x ) + £ u + y 7 

c wo 2 6 


* O* * (z w +X 0 } 

a +3 (z +x ) + Y ^ 

WO 2 
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5,6,17 Subroutine JAMES 

This subroutine provides control for the perfect gas transonic flow 
approximation computations. Newton's Method is employed to find solutions 
(s<*e subroutine NEWT) 


, sonic point displacement 
U Q > axial velocity/critical flow velocity 

a , 1st order coefficient 

6 , 2nd " 

Y , 3rd " 

to the non-linear equations (see subroutine FCALC) 



Solutions are sought such that: 

if 1st order then 3, y, f^ = o, and f^ » o are deleted 
if 2nd order then y and fg = o are deleted 

if the conical inlet angle, q., Is greater than the faring angle, then a 
special expansion is performed at the nozzle throat with u Q and T *,=o deleted. 
Expansions at points (r w ,z w )=l upstream of the faring point are performed 
with u Q and computed by subroutine FCALC and fg set zero. 


Calling Sequence 

RR r w , radial wall position coordinate 

XX z w , axial wall position coordinate 

TT Vr c » rec ^P roca ^ throat radius of curvature 

VEK k, the average gas-particle expansion coefficient 
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5.6,18 Subroutine LEGS 


Purpose 


TO solve an over de termined linear system of equations AX * b in the least 

T T 

squares sense, i.e., to solve the "normal equations" A AX » A b. 

The superscript T denotes the transposed matrix. 


Usage 

Calling Sequence: 

CALL LEGS(AB,X,BKDS,M,N,MD,ND,ATA,R,Z,SUS, SUSP, II, 12,13, 14,I5,IDLE) 


where 

AB * the M by N variably dimensioned augmented matrix (A,b) with 
maximum row dimension MD and maximum column dimension ND. 

X * the solution vector of dimension N. 

BNDS * the bounds vector of dimension N. 

M = the number of rows in the augmented matrix (A,b). 

N - the number of columns in matrix A. 

MD m the row dimension of array AB. MD i M. 

ND * the column dimension of array AB. ND £ N + 1. 

ATA = a one dimensional array in which the upper triangular part of the 
augmented normal matrix 

A T A A T b \ 
b T A b T b / 

is stored by rows. The dimension of this array mutt be greater 
or equal to (N+l) (N+"’)/2. 

R * a one dimensional array used for storage. If the inverse of a T A 
is requested, the lower triangular part is .tored by ~ows in 
this array. The dimension must be greater than cr equal to 
((N+2) (N+3)/2) - 1. 

Z c an output indicator. If the bounding -“lion is losen, L t 
indicates that the solution was affected by the be nds. 

SUS ■ the squared norm of b, i.*,, SUS « ||b|| ■ 


5 - 1 ? 7 


SUSP 


p 

■ the squared norm of AX-b, i.e., SUSP ■ | |AX-d| | 

U « Bounding option. If II • 0 the solution is unbounded and the 
routine minimizes | |AX-b 1 1 . Fbr II / 0, the routine minimizes 
||AX-b|| 2 under the side condition that 

2( X i/ B i) * 1# 

The sum extends over all i for which B A > 0. If B jL • 0, the 

routine ignores the i'th row and column of the augmented normal 
matrix and sets X i « 0. If B i < 0, the routine does not bound 

X A » The constants B ^ must be stored in ENDS. When II « 0, 

the ENDS vector is set to -1. 

12 m Solution option. If 12 ■ 0, the routine solves for X. If 

12 £ 0, the BNDS vector is placed into the X vector and SUS, 

T »1 

SUSP, and (A A) (if requested) are computed. This may be 
used for checkout purposes when the routine is used in an iterative 
process for solving nonlinear systems. .It is also used for 
comparing two "solutions" of a nearly singular system. 

13 « ATA option. If 13 « 0, the augmented norm A matrix is computed 

and placed in ATA. If 13 > 0, the normal matrix computation 
will be omitted. It is assumed the user has already calculated 
this matrix and placed it in ATA. If 13 < 0, the normal matrix 
is computed and added to whatever is in ATA. This is helpful 
in case the augmented matrix must be partitioned. 

Ik m Inverse option. If 14 / 0, the normal matrix is not inverted. 

Upon exit R will contain 

a. the lower triangular part of the matrix S, stored by rows, 
where F is an Nth order lower triangular matrix with (-1) 
on the main diagonal 

b. the solution vector X 

c. -1 

d. the main diagonal of the diagonal matrix D. 

Let C m A^A + zB ^ , then S and D have the property that SCS^* « D* 

If Ik « 0, the matrix C is inverted as C" 1 * S T D‘ 1 S and C" 1 
replaces S before exit from the routine. 

15 * ATA only option. If 15 / 0, the routine exits after forming the 
normal matrix. If 15 * 0, the routine continues by checking 12. 


* A more complete explanation is found under Functional Description. 
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IDLE * ATA row deletion option. In normal use IDLE Is a single cell 
. integer of value 0 and has no effect on the solution. The use 
of this option is best explained by an example. Suppose the 
user forms the augmented matrix and wishes to delete certain 
rows from consideration in the least squares problems (because 
of bad data or some other reason). He wishes to delete rows 
20# 31-35# 5 and l4~l6. This can be accomplished by putting the 
following sequence of integer into IDLE: 20,0,31#35#5#0,lU,l6,0 

The list is terminated by a zero in an odd position. IDLE must 
be dimensioned in the driving program if it is to be used. 


Functional Description 

1. Bounds - When the bounds option is chosen the program minim* the 
quadratic form ||AXrb|| 2 subject to the side condition £ 1, 

-2 -2 

where B is the diagonal matrix whose i’th diagonal element is B. 

if B^ > 0, whereas if B^^ < 0, the i'th diagonal element is equal to 
zero (if B^ ■ 0, the routine ignores the i'th row and column of the 
augmented normal matrix and sets ** 0). 

The method of Lagrange multipliers is applied to solve this problem 
Thus, the system of .linear equations 

(A T A + zB* 2 ) X = A T b 

is solved for various choices of the real positive multiplier z, until 
a solution X( z ) is found which satisfies the side conditions. 

The search for an appropriate value of z is systematized as follows: 

With 6 * l/lO and £g " 2/l0, the routine 

o 

a. Finds X ■ X(0). If (B~ X, X) £ 1 +6,# the problem is considered 
solved. Otherwise, 

b. define y(z) ■ (B*^ X(z), X(z)). Now, y(o) > 1 ♦ €^. Compute 
y(z, ),y(z 2 ),..., vhere z^ ■ z ^ + 10(j-l)H and H « M|X 

[B^( |A T b) i |/B ; . - A T A(i,i))] whenever this expression i6 positive, 
otherwise H « M^X [.0001 A T A(i,i)B^] 

The computation of the y(z,)*s is continued until a z is found 
such that 1 «tg £ y(zj) S 1 + ^ in which case X(z^) is accepted 
as the solution, or until two values of z^, renamed z^ and Zg 
are found such that y(z^) >1 + and y(zg) > 1 - £g. 

The required value of z is now bracketed. Now 

c. a value of z^ is chosen with Zg < z^ < z^. If 1 «£g £ y(z^) £ 1 + £. 
X(z^) is accepted as the solution. Otherwise, 
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inverse q uadr atic interpolation to zero is applied to obtain a 
new estimate z^. If 1 - < y( z^) S 1 + X(z u ) is accepted 

as the solution. Otherwise, 

e. select from the set and z h that pair which bracket the 

solution most tightly and return to (c). 

The routine will halt this iterative process if the number of solutions 
of the linear system reaches twenty. 

T -2 

2. Solution of the linear system - Let C ■ A A + zB . The routine 

T 

finds two matrices S and D such that SCS ■ D, where D is diagonal 
and S lower triangular with minus ones on the main diagonal. It 
is easy to find S and D for a 1 x 1 matrix C. Suppose they hs.ve 
been found for a k x k matrix C. Now augment C by another row and 
column 



and seek a vector CO and a scalar 3 such that 

o\ /c d^ /s T w \ /d o\ 

„ t ij (d »*; (o j ■ (o t ) ■ 


It 16 easy to verify thatuli 

B- 


S T D' 1 Sd 

2 -J 1 d 


6a + isfy the requirements. The routine builds the matrices S and D 
by the above process the final result being a decomposition of the 
augmented matrix. 

/ S 0 \ /a T A+z 3" 2 A T b W S T \ Id 0 ) 

-lj \ b T A b T b/ yo -1/ (o c</' 

The vector^? is the solution vector. The above method is a variant 
of the "square root" method. 


3* Residual Calculation - The residual is calculated as follows: 


||AX - b|| 2 - (A T AX,X) - A(A T b,X) + (b,b) 
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5.6.19 Subroutine NEWT 


Given a set of m functions in n unknowns (m * n) 

f i<V V •"* X n> i - 1, ... m 

this routine will attempt by Newton’s method to find values 

V ••• X n 

such that 


® 2 
l X 

1-1 1 


X n> 


Is a minimum. 


Restrictions 

The user must supply initial values for the solution vector and also 
a subroutine to evaluate the functions. The subroutine must communicate 
with NEWT through C0MM0N statements. Subroutine TRW LEGS is required. 

No error exits are provided. 

Method 

Newton’s method is used to iterate for a solution vector. The matrix 
inversion is performed by subroutine TRW LEGS. 

In the event Newton’s method yields a vector which is farther from a 
solution in a least squares sense than the previous estimate, the increment 
vector is halved. 


Usage 

Calling Sequence: 

CALL NEWT (M , F , MF , N , VAR , NVAR , FCALC , WF , PERTV , EPS1 ,KAXIT , 
NUMIT,BIN,AB,B,KB) 

Index to Calling Sequence: 

Input: M,MF,N, VAR, NVAR, WF, PERTV, EPS1 .MAXIT.B, KB 

Output: VAR,F,NUMIT 

Working: BIN,AB 
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Explanation of Calling Sequence: 


1) 

2 ) 

3 ) 


4 ) 

5 ) 

6 ) 


7 ) 

8 ) 


9 ) 


10 ) 


ID 

12 ) 

13 ) 


M is the total number of functions, m. 

F is the array of the m function values found by FCALC. 

MF is an array of m control words: 

if MF(I)«0, include F(I) 
if MF(I)«1, exclude F(I) 

N is the total number of variables, n. 

VAR is an array of the n variables, x^ j * l,...n 

NVAR is an array of n control words: 

if NVAR(J)*0, include VAR(J) 
if NVAR(J)*1, exclude VAR(J) 

FCALC is a name for the subroutine which calculates the 
functions, F(I). A program which calls NEWT must have an 
EXTERNAL statement containing this name. 

WF is an array of m weighting factors. These are used in 
conjunction with EPS1 to determine whether a solution has 
been reached. 

a> 1 - WF(I) 


PERTV is a perturbation factor used in calculating the 
partial derivatives required by Newton's method. If 
PERTV « s, then: 


3F 

3x 


F (x 4 6)- F (x) 
6 


_2 

where 5 = max( | ex | ,e*10 ; 


EPS1 is ac error bound used to determine whether a 
solution has been reached. If = EPS1, then a solution 
is claimed when 


m 

Z 

i*l 


W 




MAXIT is the maximum number of iterations allowed. 

NUMIT is the number of iterations required for solution. 

BIN is an array used internally by the routine and must be 
at least ot dimension n^ + 7n + 3. 
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14) AB is a tvo-dimensional array containing the augmented matrix 
of subroutine LEGS and must be at least of dimension AB (m, n+1). 

15) B is an input bounds vector of dimension n. If the bounds 
option is chosen, the iteration increments Ax are constrained 
to be less than B^ (see LEGS). B is not modified internally. 

16) KB is an input flag to be set equal to 1 if the bounds option 
is desired. Otherwise, set KB - 0. (See LEGS.) 
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5.6.20 Subroutine flNED 

The relations for one-dimensional two-phase perfect gas nozzle flow 
•re integrated through the nozzle conical inlet (Region I of Figure 5 6-2) 

The coordinate system is centered at the apex of the conical inlet. The 
differential equations are: 


dV 

-L 

dz' 


f ' r* (V - v ) 
8 p 1 

2 m * = '~ 


P r . 






J • • • j 


max 


dT 

dz* 


- 3 g ; r. 


. r 2 

P P 


J 


P C 
r P, 


(T - T ) 

_2j L 

V 

p j 


■3 " • • • j 


max 


dV 

£. 

dz 1 


V y 



V - yRT 
8 g 


v’-j 


where* : 
P 


j w 
J max n 

£ 

J"1 w 8 

8 8 


dV n 

' V V 


o m U J V A 
g g g 


* 2*r* z (1 - cos 0 i ) 2 


,2 


T ■ T 

g 2C 


V - j max w p 


-r- L ^ 


P g J’ 1 L * P J > 


C P - T ) + 


■] 


W 8 " Y 


' T iN 


V 


* The area, A is taken as a spherical cap of half angle 0 
and radius r*z'. 6 i 
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f* 

p 


J 


1+2.52 C 


1+3.78 C 


1 / 


P J P J • 


1 + 


3 . 42 C f 
J P4 


Y 1/2 P 
g * 


f (R ) (by interpolation from Table 5 • 6—*l) 
P 


2p r (V - V )/y 
g Pj g P t 

y /p r /~RT 


/ 


g g P 


j 


g 


The initial conditions are: 


200 


199 +Y 8 o 


1 + W - x > ! 


^max 

i+ E — 

■ I * 1 s_ 


W 


C 1 

P 'max p 

1 + r^S r 

V * 1 8 


C D j w 
P J max p 

1 + tl + (Y. - 1) ^ B] T -r- 1 


C 






and 


z o ■ - ir< A o /2,T(1 ' cos 9 i )> 1/2 


A 

o 


200 ' 
199+7 


i/(7-D 




The integration proceeds until: 


z' - - (1 + R (1 - cos 6.))/sin e. 

& C il 


using an integration step size 


h - Az , which is the input item DZMIN 
min 


The one-dimensional gas-particle velocity lags, 
are determined at position z* 




and thermal lags, L 


V 



The quantities 



and 

u / u* = Vg/u* 

s 

are also calculated to be used by subroutine FCALC, 


5-137 


5.6.21 Subroutine PARTH 


This subroutine is the master subprogram for the Subsonic-Trjnsonic Link. 
Tne axisymmetric gas-particle flow field of Region II, Figure 5.6-1 is computed 
using the procedure outlined below. 

STEP 1: 


The variables and are computed; 

1/2 


equil 


2C P T g 
g B 0 


/ k equil -1 \ 
k equil +1 ' 


z * -R sin 0 . -R ( 1-cos 0.) 

sci i 


R * [ 1 + R fi (1-cos 0j) ] / sin 0^ 


initially 


k 


equil 


STEP 2 

Subroutine 0NE D is used to determine conditions at the upstream 

end (z = z ) of Region II, Figure 5.6-1 
s 


Region II is next divided into zones as illustrated by the example 
given in Figure 5.6-2 JWall coordinates for the perturbation points are: 

r * 1 + R (1-cos 0^) 
w c P 

z “ R sin 0_ 
w c P 


where 


V 1 


0 


j 


— 0 . 9 ** ~ 0J» • .0, , **. 0, 

i n u i nj j 


i.e. a total of n. + n, +1 zones are constructed, 
i 3 



STEP 3: 


The perfect gas transonic flow subroutine, JAMES, is used to 
determine the variables x o , u q , a, 8, y and the transonic approximation 

coefficients for each zone. 


If faring is necessary (0^ > 0^) a special throat expansion is first 
performed and the variables u*, a*, y* and 0* are computed. 

STEP 4: 

One-dimensional gas-particle flow relations are integrated along 
the nozzle axis from z g to z ax ^ g « The differential equations integrated 

are those given in subroutine TRACE. Initial values for the gas variables 

T /T and u are evaluated using subroutine PR0P Initial values for the 
g g 0 g 

particle variables are computed from and as determined by subroutine 
0NE D. Transonic zones are exchanged as soon as 


, U Qi * U °i ± 1 

u * 2 

g 

Values K , and L. at z . are made available for use in Step 8. 
j j axis r 

STEP 5: 


A corrected estimate (constant lag) for the expansion coefficient, 
k, is computed using values at z » z a xis 


k « 1 + 


JCdL 


1 + c 2 y 


•^max p j 
° 2 " i i ~ 

j - 1 g 


w r c 


L J + K J (1 "V C 1 

L g. 


Cl w 

l + A r A l, 

\ Hi ”■ J 

1 w 
J max w oj 

1+ E — K j 
j-i * 3 
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the variable u * is recomputed as 



Using the above expansion coefficient, k, transonic conditions near the 

initial line are recomputed by subroutine JAMES. A corrected estimate 

for the nozzle mass flow, m , is found by integrating gas properties 

S 

along the initial supersonic data line using subroutine WDGI. 


STEP 6: 


Steps 2 through 5 are repeated twice to assure convergence in the 

variables k and m . 

g 


STEP 7: 

e i e i e i 

For 9 n " 20 ’ 3~ * J~ • and 6 i 
an initial position of 


r R sin 0 
o n 


z - z + R (1-cos 0 ) 
os n 


R * [l + R c (1-cos 0 i ) ] f sin 0^ 

is selected on the upstream boundary of Region II, Figure 5,6-2, and 
particle trajectories are integrated through Region II using subroutine 
TRACE. 


Tables are constructed for the following variables as found along 
the initial line: 


J » n 


® (1-cos 0 ) 

Pj n 

2 Tir*^ (1-cos 0^) 


K 


!s 


j.n u 


g 
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K 

j.n 


j,n 

2 

C j,n 


g 

T - T 
8 o P j 

"T -t 


where n - 1 , . . 4 
j - 1...J 


max 


Limiting particle streamline trajectories correspond to n ** 4. 


STEP 8: 

An initial supersonic data line is constructed consisting of 
NILP + 1 points with equal radial spacing. This data line extends from 
the axis to the wall along the parabola 


axis 


_/faxis_l Jwall j r 2 
r wall 


Limiting particle streamline points are added along this data line and 
those equally spaced points which fall too close ^within ' ^NIL P ) 


are deleted. 


Gas properties P , p , u , and v are computed at each point using 
6 g g g 


subroutine PR0P. 


Particle properties u , v , h , and p are interpolated at each 

p j p j p j p j 


point from the tables constructed in Step 4 and Step 7. The interpolation 
formulas are: 
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where 


P(x )-x + a(x. ) r 2 +b(x )r 4 +c(x )r 6 

j.n j.o J»“ j.n o,n 


c(x )- — r 
J» n iJ 1 


<r j,4 ‘ r J,3 )<X J,2 - X j,o> 


,2 x 2,2 2 w 2 2 N 

(r J,« “ r j »3 r j,2 (r j,3 " r J,2 <r j,4 " 'j^ 


(x - - x ) 
j.3 j,o 


(X J.4 " X j,o> 


2 2 2 2 2 2 
r J*3 (r J,3 " r j, 2 ^ r j»4 <*J.4 - r J, 2 > 


b(x. ) - — 5 5 

J » n ( r 2 _ r Z ) 

1 J.3 j,2 ; 


X J. Z ' X J,o 
2 


c(x. ) (r? ~ + r 2 ,) 
J.n J.Z j.-> 


a(x. „) 

J.n 


X J,2 " X j,o 


- b(x )r . 
J.n j.Z 


c(x, )r , 
J.n j.Z 


The subscript n*o denotes axis values. The variables K. and 

J t° 

L J are made available in Step 4 while K* , and \l>. are taken as: 
j»o r j.o j.o 

K' „ * k \ 1 
J.O J.l 

♦j.o-’j.l 

STEP 9: 


For each initial supersonic data line point the varables p , u , 

8 8 


v , r, z, h , p , u , and v are printed. 

8 p j p j p j p j 
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5.6.2 2 Subroutine PCALC 


The following quantities are computed: 


x ■ x , r * r , x * z 
o w w w 


dr .22 

T* " *„ <* c -*/> 


- 1/2 


d^r 2 2 9 o o 

— 9 - (R -x ) + x 2 (R 2 -x *) 

^ x 2 c w wcw 


-3/2 


.3 9 9 -3/2 7 -5/2 

“ 3x (R -x ) + 3x J (R -x l ) 

d ^3 w c w w c w 


d^r 2 2 2 2 2 4 2 2 

» 3 (R -x ; + 18x ^ (R -x 2 ) + 15x 4 (R 2 -x 2 ) 

j4 cw wcw wcw 


dx 


u* «ffx + 


u - 6 ■=■ +(b 01 +2b 00 x) r + b,, r 


2 21 22 
3 


41 


u"'* 5 V J + (2c 22 x+3c 23 x ) r + (c^c^x) r + c fil r 


u + u f + u* * + u* ' 1 
o 


2(a 20 +a 2l x) r + 4a 40 


2 ( b n x+ b 22 x 2) r + 4 O 40 +b 41 x) r 3 + bb 60 r$ 


2(c 22 x 2 +c 23 x 3 ) r + 4 (c 40 +c 41 x+c 42 x 2 ) r 3 

+ 6 (c 60 +c 61 x) r5 + 8c 80 r? 


v » + v »* + v »»» 
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and 


du 

dx 


, dr 
u + u "7“ 
x r dx 


^ - u + 2u — + u + u (j 1 ) 2 + 3u Og£-) 
.2 xx xr dx r , 2 rr \dx/ rrx ax 
dx dx 


2 

dx 2 


/dr \2 


/dr *3 . Qll dr d r 
u nr ( 3F ) + 3u rr cGT 


dx 


A 

dx 3 


u + 3u ~ + 3u — — 4 u — r 
xxx xxr dx xr , 2 r . 3 

dx dx 


dv 

dx 


. dr 

V + V “7” 

x r dx 


d 2 v 

dx 2 


A 

dx 3 


V + 2v £ + v + v f£| 

xx xr dx r dx 2 rr \ dx / 


v + 3v — + 3v + v — ~ + 3v 

xxx xxr dx xr dx 2 r dx 3 rrx 

+ v {l 1 ) 3 + 3v LJi 

rrr \dx| rr d' dx 2 


I dr \ 2 

\ dx / 


where the partial derivatives u , u , u ,u , u ,u % u . \ 

x* r xx* xr rr’ xxx* xxr* 

u , v , v . v , v , V vV .v .v , and v have been 

rrr’ x* r* xx* xr* rr' xxx* xxr* rrx* rrr 

by differentiation of the expressions for u and v given above. 


rrx’ 

obtained 
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5 . 6.23 Subroutine PRgjP 


Given the coordinate position (r,z) the gas properties u , v , 

T / T , and p are calculated from transonic coefficients previously 

g g 0 g 

determined. The velocity derivitives u , v , u , and v , may also 

8 z g z 8 r 8 r 

be computed. 

The quantity I in the calling sequence is a flag such that: 

1*1 implies u , v T /T , and p are to be calculated 
g g g g Q g 

1=2 implies u , v , u , and v , are to be calculated 

o o o a 


g z g r 


k-l 
’3 “ k+1 


q. * z_ + x - Z * z 
V 1 o w 

2 3 

P = u + a z + B + y T~ 
o o 2 6 


r = o , Properties 


u * u P 
g g o 

v =0 

f M „* 2 


g„ 


g 


P - P e =» t - 1 

g 6 ° Tg o 


r f o Properties 

* a 21 + ^21 + ^^22 + C 22 ^ z + z 

P 4 = b 41 + C A1 + 2(b 42 + C 42 ) ’ f 


P, * b, + c. + 2c. 0 z + 3c. z 
6 61 61 62 t>3 


P 8 * C 81 f 2 c 82 Z 
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P 10 “ c 101 

P 1 “ a 20 + (a 21 + b 21 ) z + (b 22 + c 22 ) ^ + c 23 z3 


a 40 + b 40 + c 40 + (b 41 + c 41 ) z + (b '-> * c ">) 2 + c,,z 

2 . 3 


P i “ b 60 + c 60 + ^ b 61 + c 61^ z + z + c ^ z 


P 7 “ b 80 + c 80 + C 81 Z + c 82 z 


62 


42 '‘W 

63 2 


2 ' r 3 
43 


P 9 “ c 100 + c 101 z 
P 11 “ c 120 
Q 4 ‘ r2 


* r 2 U 6 8 in 1 

[ p o*V *V *V * V ♦ p io' ] 
v g ■ 2 \ [ V * 2 V’ * 3 V’ * 4 v’ * 5 V’ * ‘ P ll r 

k-1 / u 2 + v 2 
g 


ii 


K-l / U TV v 


P„ “ 
g g, 


a ) 


r a o. derivatives 


G - a + Bz + y — 
O L 


u ■ u G 
g 2 g o 

v • 2 u * P, 
g r g 1 


u - 0 


v ■ 0 
g. 
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r t o , derivatives 

3 

2 


U g " “g* [ 2 ? ' 1 r + 4 P ' r3 + 6 1-5 + 8 p o r 7 + 10 P,„ r 9 


10 


V *= u 

8 , 8 . 


r 2 + 30 P c r 4 + 56 P, r 6 + 90 P„ r 8 


\-" 8 *[ 2P l* 12P 3- •'■•7 

+ 132 P 1Q r 10 1 


G 2 ■ 2(b 22 + c 22 ) + 6 C 23 Z 


G 4 = 2(b 42 + ^ + 6 C 43 Z 


G 6 m 2 c 62 + 6 C 63 2 


G 8 = 2 c 82 


6 Z g 


u g * [ G o + g 2 r2 + g 4 r4 + G e r& + V 8 ] 
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5,6.24 Subroutine TRACE 

This subroutine integrates particle trajectories through Region II 
of Figure 5.6-2 using the differential equations: 


where 


dr J 

V 

P i 


dz 

" p j 


du 

p j 

(u 8 

-v 

dz 

- “ a — 
0 


dv 

p j 

V 


dz 

- ■ a 

O u 

p j 

dT 

p j. 

-b 

0 

V 

dz 

c pt 

u 


2 V’ /J 

.N f 

A _L 

a * 
0 

2 a \ T 

P f 

; o r2 p 

L _ 

a l 


D * 
0 

a o 3 P 

r 

'■'I 


with f' and g' as defined in subroutine 0NED. 

1 J 
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Initial conditions are computed by subroutine PARTIL using subroutine 
0NED. The initial position (r^ * r^, z q ) is selected on the upstream 

boundary of Region I, Figure 5.6-1. 

The integration is terminated with the intersection of the particle 
trajectories and the end line of Region II, Figure 5. *6-1. This end line 
is the parabola 


/ Z axis Z wall \ 2 

'axis “l 2 ) r - 


wall 
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5.6.25 Subroutine WP G! 


This subroutine integrates along the initial supersonic data l'ne for 
mass flow (units of lbs mass) by trapezoidal rule: 


* “ n * 2 
* 5+1 - 4 


z i+l ‘ z i 

P u + p U - ■ (p v +0 V 1 

8 i+l 8 i+l 8 i 8 i r i+! ~ r i 8 1+1 g 1+1 g t g ± 


The Initial data line is a parabola defined as 


with 


r i+i “ r t - - 02 


z - z -I Za] tls ~ z wall 12 

i+1 axis l 2 J r i+1 

‘wall 


r l " 1 


i - 1,2,. ...50 


Particle mass flows are also computed: 


I w 


m 


\ * 






max 
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5.6.26 Program LINKS 3 


This subprogram is mainly intended to facilitate the conversion of the SPP 
code overlay structure to other computer systems. The only function of this routine 
is to call subroutine N2MAIN. 


5-153 


5.6,27 Subroutine ACOMP 


Given two adjacent points, m and n, whose properties are known, the 
Reynold's number associated with the jth particle size at a third point in 
the neighborhood is assumed to be: 


where 



A drag coefficient is found by quadratic (or linear) interpolation on 
Reynold's number from Table 5.6-1. 


f - f (Re ). 
P P 3 


A slipflow correction is applied to the drag coefficient; 

2 


(1+9. 45c) (l+2.52c)+3.03t/ 


P J p 1(1+9. 45c) (1+3. 78c)+.9U(4+11.34c)c 2 




1 + 


(r#j 


C f 


where 


(u -u + (v -v 

8 Pj 8 Pj 


Re, 


RT 


The coefficient, A , is defined as 

p j * 

J n r 

o 


f' 


P J 2 (T R) N 
8 o 


a. 


. r 2 
P P, 
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5«6«28 Subroutine ADJK 

This subroutine performs Step 4 as described in subroutine KPBPT. 
two arguments in its calling sequence, CALL ADJK(BARL , Q) are 


3ARL « X input 

Q ■ output 


The 
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5 >6. 29 Subroutine AXISPT 

This sub.?utlne uses two known points to calculate a new axis point 



Properties at point 3 are computed from properties at points 1 and 5. 
The gas and particle streamlines passing through point 3 lie along the axis. 

The following steps specify an axis point calculation. Tests are made 
to avoid particle calculations when no particles are present. 

STEP 1 and STEP 2 

Each variable at point 3 is set equal to the average of the correspond- 
ing variable at points 1 and 3 except: 



STEP 3 

Previous values of point 3 are saved, then: 
X- T A ( A r X 3 ) 



5-156 


STEP 4 


! 


U “ “ + “T* (*3 " z 5 )(B i + B i > 

P J 3 P J 5 2 3 5 J 3 J 5 


K -K + 2 < 2 3 • 2 5> (E 1 + E 1 > 

Pj 3 P J 5 2 3 5 J 3 J 5 

2p v 

P 1 P 1 
J 5 J 5 


Pj u 

J 3 Pj. 


P u 

P 1 P 1 
. J 5 J 5 


<Zo * 


* 5 > 


STEP 5 


2 <J 3 - J x ) 

V L 1 -S- G 3 


Vg ; 

P g 3 ’ J 3 “ 2 


P ■ P + 

83 g 5 


(S 3+ S 5 )(Pg 3 - P g5 )-(T 3+ T 5 )( 23 -z 5 ) 


Steps 3 through 5 are iterated until the calculated variables change 
by less than a prescribed amount, or until a maximum number of iterations 
is reached. 
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5.6.30 Subroutine CHECK 


This subroutine compares current and previous values of the point 3 


flow variables P.p.u.vh o ,, „ . 

8 8 8 8 Pj* P P j’ P d V p for relat * v e con- 

vergence by calling subroutine GRIT. If . variable's encountered failing 
the test the program returns with the convergence flag set as failed 

(N0 - 1). If convergence is achieved the subroutine returns with the 
convergence flag set as passed (N0 - 0). 
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5-6,31 Subroutine CNTRL 

The purpose of this subroutine Is to control the construction of the 
finite difference mesh for the method of characteristics solution of the 
supersonic nozzle flow. Left running characteristics are constructed 
starting at the nozzle throat point (r ■ 1, z - 0). These left running 
characteristics extend from the Initial data line or nozzle axis to the 
wall* The mesh points are calculated by subroutines INPT, AXISPT, WLPT 
and KPBPT under control of this subroutine. 

Additional points are inserted In the mesh by subroutine CNTRL by 
property averaging. Points are inserted to control the maximum mesh 
spacing. In a nozzle of shallow radius of curvature at the throat, the 
limiting particle streamlines may be very closely spaced. A point 
Insertion and drop technique allows this situation to be handled without 
the introduction of an unnecessarily large number of character 1st ice. 

No provisions have been made for the crossing of particle trajectories or 
for the occurrence of shocks. Weak shocks, however, will be ignored by 
point deletion or characteristic mapping. Figure5* 6-3fllustrates a typical 
mesh construction containing Kth particle boundary points, interior points. 
Initial line points, and points inserted and deleted. 

Whenever a mesh point calculation or insertion is completed success- 
fully, the output routine, PRINT, is called so that the point may be printed. 

Subroutine CNTRL also integrates the wall pressure by trapezoidal rule 
to determine the axial component of force existing between adjacent wall 
points, 

r 2 

4F - 2nr* 2 / P„ rdr - £(P + P )(r, 2 - r 2 )r* 2 . 

r 5 8 *■ g 2 g 5 1 5 

Total thrust is then found by 
2 12 

F. ■ F + 2irr* / P rdr 

1 i 8 

where F is the thrust across the initial line as calculated by subroutine 
C0NSTS. 
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The quantities defined below are not computed in subroutine CNTRL 
but are basic to the point calculations controlled by subroutine CNTRL. 
These quantities are given here for convenience. 


TP. 


- 


/p (u 2 + v 2 ) - Y? 
g g g g 




u p - yP 
g g g 

p u v + .R -R 

&._B— 8 L- 2 - 


p g u g v g ~ 1 R 2 R 
1 V 
V 

-R 

u 

g 


P< u 


P4 

j 

p 'jN (u - u ) 

!& _E Ii_ 

P u 

■ gj p i 

f (v - v ) 

>!*! 1 !i. 


1} 


j 


it 

P 8J 


yP V 
t,t - — 


° g (v g 1 R + U g 2 R) 

1 R x v 


°* (V » 1» ~ U * 2R) 
X R 1 V 
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*3 “ *2 1 

Pg + -- - 2 - (P 3 u 3 + F 2 y 2 ) + f( 3 R 2 + 3 R 3 ) 


u 


(6 8 

{ ■ ~r < p g 2 v g 2 + p g 3 v g 3 > + it k 2 u g 2 + p g 3 u g 3 > 


(r 3 - r 2 ) 


(* 3 " z 2 ) 


p 

N 

P 

8 J 

f + 

g 2 

p , 

2 3 

P 

g 3 


«2 


N 


2*2 J 


P 

®3 

N 

P 

g 2 

i P J» J 

g 3 

3*3 + 

k J 

8 2 


N 


3 t 2 l 


+ 4 t« g2 (r 3 - r 2 ) - y^z 3 - zj 


\ '2 * *., *2> 
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'p ^ 
8 


i. p » J 

g 3 


l f 3 + 


lV } 


<z 3 " z i ) 1 

P gi < F 3 P 3 + F 1 V - I<3 R 1 + 3 V 




V + p„ v ) + -T— (p u + p u ) 
8 1 g 3 8 3 2 8 1 8 1 8 3 8 3 


(r 3 - r l> 


(Z 3 - 2 1 ) 


P 

8 3 

N 

P 

8 1 

l P o J 

8 3 

2*3 + 

k 

h 

P 

g 3 

N 

f 

P 

8 1 

k J 

g 3 

3*3 + 

l° 8 i 

g l 


2 e ll 


N 


3 t l J 


+ T Iu 8l (r 3 - V - v 8l (z 3 - z l> - u g 3 r l + v g 3 z ll 


fP 1 

8 ’ 

N 

+ 

N 

•1 

11 

° g 3^ 
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Note: When calculating an axis point the evaluation of reoulres an 

Indeterminate (v - r * 0) quantity, . In this case the 
S ^ 

expression below is used: 


► « 

V 

u 

V 

i 

00 

g 3 

r 

* a 

3 r l 

V 

8 1 

L u + (z, - z,) - — J 




v 

This calculation amounts to extrapolating for the ratio, assuming 
that the flow near the axis is a source flow. 
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5.6.32 Subroutine CONSTS 


An Integration Is performed using trapezoidal rule for thrust 
along the initial line: 


F = 2„r*2 


j- r 

-wall ( -v ) + T, 0 u (u -v ) rdr, 

J [ g H g g g r j i p Pj Pj Pj Pj 


where K is the total number of particle sizes present at the point under 
consideration. 

The indeterminate quantity 

v 

r 

at the initial line axis point, m, is calculated as 


& 


lx ) 

r j 


^m-1 

Vl 


g, 


m 


L Vl 


+ (z_ - z_ ,) 


°m-l 


m m-1 r 


m-1 J 


This quantity Is required in both the axis and off axis point cal- 
culations. 
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5,6.33 Subroutine CRIT 


The purpose of this subroutine is to compare tvo values for relative 
convergence and return a flag stating if convergence has been achieved. 

Calling Sequence : 


XN 

m 

NO 


m 


A flag such that 


if x - 0 then NO - 0 
n 

if |(x_ - x )/x |s 5.10’ 5 then NO * 0 
m n n 

if |(x - x )/x |> 5.10“ 5 then NO « 1 

m n n 
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5.6.34 Subroutine CUBIC (X. Y. VP, N, ARG, YARG) 

The purpose of this subroutine is to perform cubic interpolation 
for a tabulate 1 function whose derivatives are known. 

Calling Sequence; 

X 

Y 

YP 

.N 

ARG 
YARG 

Restrictions : 

The calling program must define arrays for the dummy variables 
X, Y # and YP. These arrays must be at least of lengths N + i 9 N # and 
N respectively. The subroutine will save its last used table position 
number in X(N + 1). 

If x < the program returns y = y^ 

If x > Xjj the program returns y = y^ 

Mfilhgd: 

Given: x ^x<x, 
o 1 

so that y Q# y o ' # y^, and y^ 1 are known. The cubic interpolation formula 
given below is used to determine y 

y = A(x - + B(x • X Q )^ + C(x - x q ) + D 

where 

a = ^ ((y i' + y 0 ‘) h - 2k I 

B -prlfri +2 ^. , ) h - 3k I 


independent variable 9 x^ 9 such that 

dependent variable 9 y., = y(x^} 
derivitives of the dependent variable 9 

Yi s 

is the number of entries in each of the above tables; 
is i # ... N 

is the argufttent, x 9 for which interpolation is requested 
is the result 9 y = y(x) 


is a table of the 

*i+l ix i 

is a table of the 

is a table of the 
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D s y 

• J o 

ka U-r 0 
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5.6-35 Function EF N 


The quantity is calculated In this function subroutine. A test is 
made on particle enthalpy, h , in order to determine the proper particle 


phase. 




v 


u 


N 




u 


N 


K * V 1 ■ " 

I 1 -. - r%] 


• b p R) g 1 r 
p i I p / * 

J \ at 


h - h p 

Is Ll_ A. 

Cp R p 


h > h 

p i p / 


h 4 \ < h 
p 8 Pj 


h < h 


Pj Ps 
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5,6,38 Subroutine ERR0R 


The purpose of this subroutine is to print error messages as requested 
by subroutines COTRL, ADJK, and WLPT. See Appendix A for conversion aides. 
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5.6.37 Subroutine KP3PT 


This subroutine uses five known points to calculate a Kth particle 
boundary. 



Properties at point 3 are computed from properties at point 2 and 
interpolated properties at point 1. 

The following steps specify a Kth particle boundary point calculation. 
Tests are made to simplify the calculations when a gas and one r article 
boundary point is to be computed. 

Step 1 

Each variable at point 3 is set equal to the average of the corres- 
ponding variable at points 2 and 4. Point 1 is initially set equal to 
point 4. 
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Then 



J ■ If •••* 
j ■ If • • • f K-l 


Step 2 

Previous velues of point 3 are saved, then 


7 - ta(< 2> *3) 



TA(C , C ) 
\ K 3 



r 3 “ r z k ( z 3 " z 2^ 


Step 3 

I - ta(x 3> x x ) 

Points 6 and 7 are selected such that point 6 is point 4 and point 7 
is the next point on the previous characteristic. 
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Step 4 


This step is performed by subroutine ADJK. 




If q x > 1 
If q x < 0 
If 1 > q x > 0 


select new points 6 and 7 which are one point 
higher on the previous left running characteristic 
and restart Step 4. 

select new points 6 and 7 which are one point lower 
on the previous left running characteristic and 
restart Step 4. 
proceed . 


Step 5 




r 6 + 1l< r 7 - r 6 } 


P + q. (P - P ) 

8 6 1 8 7 8 6 

P o + - P » > 

g 6 1 g 7 8 6 

u + q- (u - u ) 

8 6 1 8 7 g 6 


v + q, (v - v ) 
g 6 1 s 7 g 6 


h + q,(h - h ) 
P< 1 P< Pj 




% 


V + - V > 

J 6 J 7 J6 


j ■ 1 K-l 

J “ 1 K-l 


u + q, (u - u ) 
P 4 P 4 P 4 


j • 1, .... K-l 


v - v + q, (>r - v ) 

'll X 1 X X 


j * l t • • • i K-l 


Step 6 

Steps 6 and 7 are performed for each j - 1, . . . , K. The gas prop- 
erties at point 4 are first saved. 


If J - K 


If j * K 


restore gas properties at point 4 and go to 

Step 7. 

proceed. 


TA(C , C ) 

p 4 p 4 

h h 


r 3 + V (Z 2 " Z 3 - r 2 T, 


1 - C 


, Z 1~ Z 2 
' P j r l " r 2 


r 4 ~ r 2 
r l “ r 2 


Z 4 “ z 2 + q 4 (z l " z 2 } 


" + q 4 (P g " P « > 

g 2 g l g 2 

- ♦ q 4 (» gi - », 2 > 

■ \ * - v 

• \ * v« 4l - y 


V * "‘'v - V > 

J 2 3 1 J 2 
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Step 7 


Step 8 
for j - K 


V • V + ’4<V - v > 

3 4 3 2 3 1 3 2 


u - u + q. (u - u ) 
P J4 P J, 4 P 3i P J, 


■ v 4* q. (v v ) 

p i 4 P 1 p 1 
3 4 3 2 3 1 J 2 


•it "j 


u ?< + iS- v<vV 


V + 2 1 ^3- V\ + \) 

J 4 


V *7"3 -V<VV 

3 4 


3 J 4 

+ E ) 
3 3 4 


Pj 3 l V <r 3 2 “ r 2 2) " v p ( 2 t- z 3 ) (r 2 


r j- 


r J. 


I P [u (r 2 - r 2 ) 

+ r 7 )l ' P J 4 P J 4 4 2 


- V ^-V (r 2 + VI + V [U P (r « - r 3 > 

3 4 J 2 ^2 

- V I( *4 - *2 )(r 4 + r 2> - <*3 - *2 )(r 2 + r 3 )1] } 

J o 


for J < K 


V (r l 2 - r 2 2) 

3 3 


J — { 

- V [U 3 - Z2 )(r 2 + r 3 ) + (z x - z,)^ + r 3 >] { 

J 1 


P n [u (r 

P 1 P 1 
J4 J4 


(r 3 - ) - v [(z 4 - z 2 )(r 4 + + 0^ - z 4 >(r 4 + r )]] 
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+ \ lU P 4 ^ “ r 4 2) “ V I(Z 1 ' Z 4 )(r 4 + r l> " (Z 1 - z 3 )(r l + r 3 )]1 
3 1 3 1 

+ P p 4 f V (r 4 2 ■ r 3 2) - V t(z 4 - z 2 )(r 4 + V - (z 3 “ z 2 )(r 3 + r 2 )]1 } 
3 2 3 2 3 2 

Step 9 


C 

g 


r 


5 


TA(C , C ) 
g 3 g 5 


r 3 + C g [Z 2 - Z 3 ' r 2 71 


- Z 1 " Z 2 
1 - C — - 




r 5 - r 2 
r l " r 2 


Z 2 + q 5 (z l ‘ z 2> 


P + q, (P - P ) 
P g 2 + q S (p 8l - 


u + q c (u - u ) 
g 2 5 g i g 2 


V K + (v g ' V b ) 
®2 ' 8 1 g 2 


h + q,(h - h ) 

P J 2 5 % P J 2 

P p + q 5 (p p ’ c p ) 

J 2 J 2 

U + q e (u - u ) 

P-i j P4 Pj 

J 2 J 2 


v + q c (v - v ) 

p j, 5 p j, p j. 


j “ 1 K-l 

J - 1 K-l 

J “ 1 X-l 

.1 - 1 K-l 
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A test Is cade to determine whether point 5 Is above or below the 
Xth particle streamline. 

If 



then 

h - p m u *= 0 and primes must be added in 

% % *5 


Step 

10 

to J v G 3 , Gj, Hj, Hj, T 3 , and Tjj 

otherwise: 

h 


y(h 

+ h ) 

% 


2 P K 

4 

P K 

*2 

P 

a. 

f(p 

+ P ) 

% 


2 P K 
*4 

P K 

*2 

u 

m 

“(n 

+ u ) 

P Kc 


2 P K, 

P K 0 

5 


4 

2 

V 

m 

|(v 

+ v ) 

% 


% 

P K 

K 2 


Step 10 

For the primed quantities the implied summations are to be taken for 
j ■ 1, . . . , K-l, not K. 

2(Jj - J 2 MQ 3 + Q 2 + h 5 + a 3 ) - 2(J 1 - j 2 hq* + q 3 + Qj + q 2 ) 

U g 3 “ (L’+ L 3 + L’+ L 2 )(Q 3 + Q,+ Hj+ H 3 )-(L 3 + L 2 + G 5 + G 3 )(Q’+ Q 3 + Qj+ Q 2 > 



2(j; - J 2 ) - u g (Lj + L 3 + + L 2 ) 

Qj + Q 3 + + Q 2 


J 2 + 


U g 3 (L 3 + V + + V 


2 


<S 3 + 


Sq 


- P 


S 5 )(P g 3 


" V " 


(T, 


+ T 5 )(*3 - 


V 


8 « 


Steps 2 thr ou gh 10 are Iterated until the calculated variables change 
by less than a prescribed amount or until a maximum number of Iterations is 
reached. 



5,6,38 Subroutine N2MAIN 


The purpose of this subroutine is to call the two portions of the axisymmetric 
supersonic method of characteristics subprogram (subroutines N3MAIN and CNTRL) . 


5,6,39 Subroutine N3MAIN 

The purpose of this subroutine is to call subroutines C0NSTS and WALL. 
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5.6,40 Subroutine PRINT 


The purpose of this subroutine Is to print output for the supersonic 
flow calculation as requested by the program input. Printout may occur 
after the completion of each mesh point calculation. Points for print are 
selected as follows: 

1) The following points are always printed: 

axis points 

Kth particle boundary poluts 
initial line points 
wall points. 

2) Interior points are selected for print only along every n^th 
left running characteristic and only at every n^th position 
along these characteristics. 

3) Inserted points are printed if all points are to be 
printed (n^ * n 0 ■ 1). 

The items printed are listed in the table below in the order they 
appear, left to right, on the output sheet. A header is printed for 
identification purposes above each characteristic. 

Row One: 


Item 

Header 

Meaning 

Units 

L.R.C. number 

LRC 

left running characteristic number 

none 

Xdent. number 

ID 

type of point (see below) 

none 

r 

R 

r position coordinate 

none 

z 

Z 

z position coordinate 

none 

M 

MACH 

Mach number 

none 

T 

g 

TG 

gas temperature 

°R 

V 

ft 

VG 

gas velocity (scalar) 

ft/sec 
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Row One: (Cont'd) 


Item 

Header 

Meaning 

Units 

6 

g 

THETA-G 

streamline angle 

degrees 

T /T 

8 *0 

TG/TGO 

ratio of gas temperature to chamber 
temperature 

none 

p /p 

8 8 0 

PG/PGO 

ratio of gas pressure to chamber 
pressure 

none 

P„/P 0 

8 8 o 

DG/DGO 

ratio of gas density to chamber 
density 

none 

V* P /P 

L ?J 8 

SDK/DG 

ratio of total particle density to 
gas density 

none 

C F 

CF 

thrust coefficient 

none 

I 

sp 

ISP 

specific impulse 

sec. 

iteration 10 . 

IT 

number of iterations required 

none 


Rows Two through K+l 

A row is printed for each particle size, k>l,...,K. 


Item 

Header 

Meaning 

Units 

k 

K 

particle size number 

none 

**k 

REK 

particle Reynolds number 

none 

V 

Pi, 

VPK 

particle velocity (scalar) 

ft/sec 

K 

e 

Pk 

THETA-K 

particle streamline angle 

degrees 

K 

T 

P k 

P„ /p„ 

P k 8 

TPK 

particle temperature 

°R 

DPK/DG 

ratio of particle density to gas 
density 

none 

P n / 0 r, 

p k p o 

DPK/DPO 

ratio of particle density to chamber 
particle density 

none 

r 

Pv 

RPK 

particle radius 

ft. 


■i 
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Of the quantities above, the following are calculated within sub- 


routine PRINT: 


T , " p g R 


V * * u + v 
g g g 


e - ate tan | • 57.29578 

8 \V 


2 P r . /p T R 
« P k 8 8 0 

■■k* ”57- 

8 0 \ 8 


(u - u ) 2 + (v - V ) 2 

8 Pfc 8 p k 


k-1,.. .K 


'* *2 

tr p. 


SP w + * v 

8 E Pv 


h - h 

T - T + 

p k p m C P / 


for h > h 
p k p £ 


T - T 
p k p m 


, for h ^ h ^ h 
p s p k p / 


p k C P 


, for h < h 

p k p s 


5-182 


The table below relates the printed Identification number to the type 


of point calculated: 

ID 

1 

2 

3 

4 

5 

6 
7 


Type of Point Calculated 

initial line point 
interior point 
axis point 

Kth particle boundary point 
wall point 

point inserted on the previous L.R.C, 
point inserted on the R.R*C. 


In addition, whenever a wall point is printed, subroutine ISP1D is called to 
calculate the ID reference ISP and the quantity, is calculated as 


I (r ,Z) 

sp T D2P w 
*TD2P “ I (r^J 

s PlDO A w 


5.6.41 Subroutine PUN T 


This subroutine uses three known points to calculate an Interior 
point. 



Properties at point 3 are computed from properties at points 1, 
2, and 4. Point 4' is required only for the calculation of particle 
density, P pj^» ** point 4' is not available (IFG ■ 0) an alternate 
calculation is used (see STEP 6). 


The following steps specify an interior point calculation. Tests 
are made to avoid particle calculations when no particles are present. 

STEP 1 

Each variable at point 3 is set equal to the average of the 
corresponding variable at points 1 and 2. Other values required for the 
first iterations are set as: 




3 

3 


J“l> • • • »K 
j-1 K 
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* 


STEP 2 


Previous values of point 3 are saved, then 
k • T A ( < 2 $ 

7 - T A ( A x , X 3 ) 


r l - r 2 + < *2 “ * z i 


< - X 


STEP 3 


t 3 m r 2 + k (z 3 - z 2 ) 


c r IA <v v 




«r*2 


z 2 " z 3 ' r 2 r 1 -r 2 


r 5 ‘ 


V*2 


1- C 


8 r r r 2 


r 5 -r 2 


r l _r 2 


*5 “ z 2 + <>5 ( V Z 2 ) 

P - P -0 q. (P - P ) 
g 5 8 2 5 8 X 8 2 


P g 5 “ P 8 2 + q 5 ( P 8l - P g 2 ) 
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\ ■ \ * 'i\ ■ v 


V o ’ V „ + 1s< v . " ▼_ ) 

85 82 ^ 8 ^ 82 


h - h + q (h - h ) 
P« P 1 5 P 1 p < 
•>5 J2 Jl J 2 


J-l.. 


P- * P„ + «Jc(P„ - P- ' 

P< P< 5 p. p 4 


’5 • 2 


h % 


J-l.. 


VI ■ U + q c (u - u ) 

p i p i 5 p i p< 

j 5 j 2 Ji j 2 


j-i.. 


V n " v „ + < U< V „ - v - ) J-l... 


P J P 1 5 p 1 p 1 
J 5 J 2 Jl J 2 


STEP 4 


Steps 4, 5, and 6 are repeated for each j=l K. 


C - T A (C , C 

P J P 1 P 1 

J J 3 I 4 


r, + C 


T 2 1 ~ z 2 1 

3 ' ^jL' 2 ' 23 " f2 *i - r 2 J 


- Z 1 ~ Z 2 

1 - C — 1 


P J r l ' r 2 


4 r l " r 2 


*4 ’ *2 + q 4 (z l " z 2 } 


P„ - P + q, (P - P ) 

84 g 2 « g x g 2 
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\ - X * "*'% - V 


u • U + fl. (u - u ) 
«4 *2 4 *1 *2 


V “V +Q/(v - V ) 

*4 *2 4 *1 g 2 


h p 4 “ V + q 4 (h P . ■ h P , ) iml K 

3 4 3 2 3 1 3 2 


p Pl " p p, + vv • p Pl > J-1 K 

3 4 J 2 Jl J2 


u -u + q.(u “ u ) j-l,...,K 

P* p i 4 P 1 p i 

3 4 J 2 3 1 J 2 


v -v + q. (v -v ) 1-1,..., K 

P, P. H A V p p 

h J 2 h 


STEP 3 


Calculate 


V - V + — <*3- V <V +B p > 

J 3 h j 3 J 4 

A 

P J 

v - v + -r (*« - z )(D + D ) 

P 1 P 1 2 3 A Pi 

J 3 J 4 J 3 3 4 


V “V + 2 ( * 3 -V ( V + E p > 


’3 J 4 


h h 


5-187 


STEP 6 


If point 4 ' is not available (IFG • 0) : 


,, "J - r 2> ' V [ r 3 - r 2 * "l + ' 

J 3 3 3 L S 


3 )(z r z 3 ] 


\ V (r l ■ r 3 )_ V I (Z 1 " *2 )(r l + r 2> 
J 2 J2 J 2 L 


2 2 
r 3 " r 2 


/ 2 2 v 

V r 3 - r 2> 


-v p (z r z 2 )(r 1+ r 2 ) - (r 1+ r 2 )(z r z 3 ) 
Ji L 


If point 4 * is available (IFG>0): 


**^3 u (r l - A 
% 1 2 


[ (* 3 -* 2 ) (r 2 +r 3 ) + (z r z 3 )(r 1+ r 3 )] 

Jo 


u Pj (r r r 2 )-v Pj [ ( V' Z 2 )( V +r 2 ) + < z r*4 ,)(r 4 ,+r l ) " 


> u ( r i -r 2 >- v D 

P j 4 . P J 4 . P J 4 . 


L 

.. V (r H')- v p 4 [ ( V , 4' )(r *' +r l ) - < *l-*3 )(r l +r 3 ) ] 
J 1 J 1 3 1 

- 1 


+P_ u (r.,-r?)-v 
P 1 P 1 ^ 3 P 1 

J 2 J 2 J; 


[ < z 4 <" z 2 )(r 4 ,+r 2 ) " (z 3‘ z 2 )(r 3 +r 2 ) ] 


STEP 7 


2(J 3 -J 2 ) (Q 3 -H) 2 +H 5 H 3 )-2 (2Q 3 -H3 1 +0 2 ) 


5 3 tfLj+L^) (Q 3 +<3 2 +H 5 +H 3 )-(L 3 +L 2 +G 5 +G ? ) (2Q 3 +Q 1 +Q 2 ) 
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2(J 3 -J 2 )-u g (2L 3 +L x +L 2 ) 


V 8 3 “ 2Q 3 ^ 2 




P„ " P„ + 
*3 8 5 


“g 3 

2 

( W ( ? , i -F, < )-(i J « 5 )(, 3 ., 5 > 

2 


Step 2 through Step 7 are iterated until the calculated variables 
change by less than a prescribed amount, or until a maximum number of 
iterations is reached. 
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5.6.42 Subroutine SUMP1 


Routine SUMP1 calculates the qualities i t where 


1* » T) A P 

j-l P J P J 


2 ‘ - E v v u . 


j=i p j p j ? j 


f' A P v 

j=i P J P J P J 


k r j 2 

(y-i) z A n p „ (U ■ -U ) + (v -V ) - 

J-1 p j P J L 8 v i J 


C E u 
p 

3 Pr 


u 1 
i-il 

(if 

iW 


5.6.43 Subroutine SUMP2 

Subroutine SUMP2 calculates the same quantities as described above in 
SUMP1 . Routine SUMP2 also calculates the additional quantities l* where K is 
replaced by K-l, which are required by the particle boundary point calculation. 
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5-6.44 Subroutine TAFN 


This subroutine calculates an angle jveraged tangent using the 
relation: 


r* 


TA ■ tan I *r (arc tan X + arc tan 


Y) ] 
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5,6,45 Subroutine WALL 


This subroutine makes available five options for the specification of nozzle 
wall contours. Either a cone, parabola, arc, contour, or, a spline fit curve may 
be attached tangentially to a circular arc. The circular arc is specified as follows: 

Symbol Definition 

Angle defining the downstream end of the arc 
0j Angle defining the upstream end of the arc 

r/r* Normalized arc radius 

where 



One hundred points are generated along the arc for equally incremented 
values of angle, q: 

r = 1 + r/r* (1 - cos q) , 0 L = 9* 0j 
x = r/r* sin ^ 

Five options are specified below for attaching wall contours tangetially 
to this circular arc. The cone option requires the nozzle expansion ratio as In- 
put. The parabola and arc options require the exit point coordinates, r and x , 

v C 

as input. The spline fit option requires a smooth table of wall coordinates, r { 
and x., and the exit angle at the nozzle lip, o , as input. 

j e 
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Option 1, Cone 




x c = (r e - r a )ctn 0 + x g 

Option 2, Parabola 

Two hundred points are generated at equally incremented values of x such 

that: 

r = a + v b(x - c) 


where 


*% - 1 ~ 2r a (x e ~ x a } tan 
2 C r e ' r a ' U e " V tan »j] 


b = 2 (r - a) tan 


c = x - 
e 


(r - a) 2 
e 


Option 3, Arc 

Two hundred points are generated at equally incrmented values of angle 
such that: 

r = r a + R~(cos q - cos e^) 
x = x a + R(sin 0 ^ - sin q) 

where 

R = 


(x - x ) 2 + (x - r ) 2 
e a e a 


T[(x e - x a ) sin 6j - (r e - r a ) cos ^ 
and the increment on angle is 


200 


where 


»e = Sin 


■*[- 


sin 9j - 


X (x e - x a J 


'] 


Option 4 , Spline Fit 

Two hundred points will be placed along a curve input in tabular ft mi by 
a cubic spline fit method (see subroutines CUBIC and SLP) . Derivatives are cal- 
culated from the input values for 0. and the exit nozzle lip angle, 0 . This curve 

J 6 

must connect tangentially to the circular arc described above. 
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5.6.46 Subroutine WLPT 


This subroutine uses two known points to calculate a wall point. 



Properties at point 3 are computed from properties at points 2 and 
5. The wall is necessarily a gas streamline. It is assumed that no 
particles are present at points 2 and 5. 

STEP 1 

Each variable at point 3 is set equal to the average of the 
corresponding variable at points 2 and 5. 

STEP 2 

Previous values of point 3 are saved, then 
< * T A (< 2> »c 3 ). 


STEP 3 
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If 


«3 > 1 + C w 
If 

^3 < " e w 


If 


select new points 6 and 7 upstream 
one point and restart STEP 3. 

select new points 6 and 7 downstream 
one point and restart STEP 3. 


- e w < fl 3 £ l + e w proceed. 

r 3 - r fi + q 3 (r ? - r 6 > 


STEP 4 



2(J X - J 2 ) 

L 3 + L 2 + Gj + G 3 + C g ^ (Q 3 + Q 2 + » 5 + H 3 ) 


v ■ u C 
8 3 8 3 8 3 


P - J, - 
8 * 1 


u (G. + G,) + v (H, + H,) 
g 3 5 3 g 3 5 3 

~ - 

(S, + S c ) {? - ) 


go g 


+ 


3 y w g 3 g 5 


Steps 2 through 4 are iterated until v e calculated variables change 
by less than a prescribed amount, or until a uia::imum number of iterations 
is reached. 
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5.6.47 Subroutine XS1.P 

This subroutine Is identical to subroutine SLP (Section 5.5.15). 


5*7 Link 60 Subroutines 


This overlay and subroutines contain the TBL Module. The subroutine 
write-ups are from Reference 5 with few exceptions. 

Only a minimum of modifications have been made to the TBL module. These 
modifications deal mainly with the linkages between the various modules of the 
SPP code, 0DE and TD2P, and the TBL module. 


5.7.1 Program Link 60 

This subprogram is mainly intended to facilitate the conversion of the 
SPP code overlay structure to other computer systems. The only function of this 
routine is to call subroutine TBL. 
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5.7,2 Subroutine BARCON 


Subroutine BARCON accepts initial conditions and executes a Runge-Kutta 

solution for the boundary layer along the contour. The initial conditions are set: 

0= 0 Q , e = 9 0 , at x=x 0 . A step-size Ax, used to increment the axial coordinate x, 

is determined from input quantity Ax w=v and the local entires in the x-table. Having 

in ax i i 

two first order ordinary differential equations dg/dx = f(x, 0, 0) and d0/dx = g(x,9,<fl) 
and using subroutine BARPRO to evaluate the derivatives, a four term Runge-Kutta 
numerical solution is used. 

1. Evaluate 

f i = _ a£l ' 9 i = !xl 

x^ x„ 

o o 

Set 

9 = 9 o + T" (f l J ' $ = + fep 

2. Evaluate 

f 2 = ^ Xo+ Ax. ' g 2“dFl Xo + ^c 
Set 

Q=0 o + “^r^V # + T” k 2 ) 


3. Evaluate 


Set 



0- + Ax (f 3 ) , 0 = 0 O + Ax(g 3 ) 


4. Evaluate 


f 


_ d0 | 
4 "3x ' 


x. 


+ ax 



+ Ax 


5. At x + ax, 0 and 0 are then evaluated using 

0 = 9 o +i f (f 1 + 2f 2 + 2f 3 + f 4 J 
0 > = *0 + (g l + 2g 2 + 2g 3 = 9 4 ) 
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Now using these values of x, $ , 0 as the Initial values for the next Interval, we 
repeat the above procedure on through to the end of the defined contour. 

This subroutine also calculates and prints a table of normalized contour 
points corrected for displacement thickness. This print Is made to assist in the 
preparation of a potential flow nozzle contour table for input ot the TD2P, TDK, 
or other method of characteristic programs. This dimensionless potential flow 
(adjusted) nozzle contour Is calculated by the following relationships. 

POTENTIAL = y t * y t,REAL " 5 t 


^POTENTIAL = [y * y t(REAL - 6 * cos or]/y^ POTENTIAL - 

Potential = ' y t, real + 6 * sln « 3/y tj potential 

where 


y t,REAL - actual (geometrical) throat radius 
x are XI TAB 
y are YITAB 


oi= arc tan dy/dx. 

The nozzle throat is indicated by the subscript t. 
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5.7.3 Subroutine BARPRO 


Subroutine BARPRO calculates the boundary layer, given energy thickness 
( 0 ) and momentum thickness ( 0 ) and d 0 /dx and d 0 /dx along the contour at a point 
(x). The interpolation routine XNTERP is used to evaluate inviscid flow and contour 
properties and derivatives at point x. Subroutine ZETAIT is then entered, and the 
parameter f and boundary layer thickness A/ 6 # 6 * are computed. An iteration pro- 
cedure is then used to calculate skin friction coefficient (C^) and Stanton number 

<°B>- 


1 . 

2 . 


An initial guess (C, a ) is made 

f g 

T 

Calculate C f R- = (^pL) X " m C fa R 0 


3. 

4. 


5. 


Subroutine CFEVAL is used to evaluate C f as a function of C^R^- 


Calculate (-y— ) and C fQ = — 


'i 


__ T m 

aw ^ aw j / s j 


— j — ) 

e aw 


C fa- C fa 


A relative error comparison Is made. If | — ^ — — 2— | < TOLCFA 

fa g 

convergence is considered to be attained. If unconverged, a form of 
Wegstein's method is used to calculate a new guess (C^ ) and steps 2 
through 5 are repeated up to a maximum of 50 iterations. 9 


The derivatives d 0 /dx and d 0 /dx are now evaluated, from the differential 
equations, at the point x. For the point x at the end of a Runge-Kutta step (IND=1), 
total heat transfer and drag quantities are also computed. 


This subroutine also prints the program output under control of the input 
variable, IPRINT. 

This subroutine was also modified to print and write (punch) out the 
linkage data needed by the SPP master control module. In particular this routine 
calculates the force decrement due to a turbulent boundary layer using the equation 


AF 


(2iry) k o e 



0 • COS Cl (1 


e 



) 


-1 dv 

where: <* = tan -gj- 

K = 0 for planar flow 
K = 1 for axlsymetric flow 

and the decrement in I is calculated from 

sp 

^sp ' aF/l " 

where: m = mass flow from input or the TD2P module 

= force decrement as calculated above 
The above items are only printed out: when the wall slope is positive. 
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5,7,4 Subroutine BARSET 


This subroutine coinputes program constants and sets up C and H versus 

P 

T tables and tables of P and T versus x for the inviscld flow. Based on the input 
table of C versus T, subroutine BMTAB is used to generate coefficients of a 

r 

cubic spline fit. These coefficients are then used to compute piecewise partial 
integral values to be used by subroutine SEVAL which relates C^, T, H, S and P. 
Subroutine GETPT is then used to evaluate pressure and temperature at each point 
of the Mach number table. 


5.7.5 Subroutine BMTAB 

This subroutine computes the coefficients of a piecewise cubic spline fit 
of tabular data , 


5.7.6 Subroutine BMFIT 

This subroutine computes the coefficients of a piecewise cubic spline fit 
of tabular data. 
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5,7.7 Subroutine CFEVAL 


Subroutine CFEVAL evaluates the adiabatic skin friction coefficient, C^, 
as a function of C^R-, based on Coles' relationship. For values of C^R— , such 
that 0 <C f Rg <2.51, C f is evaluated from 

~ .009896 

f (c f v * 562 

For values of C^R— , such that 2.51 < C^R— < 1982759.2, C^ is evaluated from a 
piecewise cubic spline fit of log (C^) vs log(C f R^). All other values of C^ 
are considered out of range. 


5.7.8 Subroutine DIRECT 

Subroutine DIRECT controls the overall flow of the program. It calls sub- 
routines which successively read in data, compute program constants, fit curves 
and then solve the boundary layer solution algorithm. 


5.7.9 Subroutine EDITCP 

This subroutine attempts to edit a table of C p versus T from the boundary layer 
edge condition table of C p versus T if the IFEDGE=1 or 2 option was selected in 
the $TBL namelist and if n<? C versus T table was input. 

r 
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5 f 7 . 10 Subroutine FIIF 

Subroutine FIFF defines the function (of "s") to evaluated by the numerical 
integration routine INTZET. Three second degree coefficients (Ap._Bp, C p ) and an 
exponent value (m f ) are set by subroutine ZETAIT. Enthalpy value h is computed 

from 

h = Ap + Bps + CpS® 

Entering subroutine SEVAL with h. we obtain!. The appropriate form of the function 
is then evaluated: 

For I P = 1, f(s) = (-r) s m f (1 - s) 

F t 

I„ - 2. f(s) - (il s »« 
f t 

where T/r has been used for q/q . 
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5 .7.11 Subroutine GETPT 


For a given Mach number, subroutine GETPT computes pressure P and 
temperature T e from . H and T relations. For constant specific heat (y = y Q ) 



i y “ 1 

(1 


For variable C p the following iteration is used 

1. An initial guess is made for temperature T from 

g 


T 

g 



M 2 


2. Using subroutine SEVAL, C p and H are computed at temperature T 

3. The specific heat ratio y is then evaluated from 


C 



4. Temperature T is then calculated 

2(H - H) 

o 

e y RM 2 

5. A relative error comparison is then made. If 


T - T 

!-^P — 2-' < TOLZME 

g 

convergence is considered to be attainea and subroutine SEVAL is used to evaluate 
pressure P 0 at temperature T g . If the convergence criterion is not satisfied, a 
form of Wegstein's method is used to calculate a new guess for T g and steps 2 
through 5 are repeated up to a maximum count of 50 iterations. 
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5.7.12 Subroutine INTZET 


Subroutine INTZET, given a function definition and limits of Integration 
(x^ to Xg), uses a numerical technique based on the qulntlc spline to evaluate the 
Integral of the function from to x 2 . The Interval over which the Integration is 
to be performed Is divided Into 20 coarse subintervals. Since the limits are re- 
stricted to certain values (0, 1, £ , l/£ ), a check on the length of Interval (x 2 - 
Xj) Is made. If (x 2 - Xj) < .10, the function Is evaluated at the 21 sublnter/al 
endpoints and the Integral Is approximated by the qulntlc technique, with a "one- 
sided Simpson's Rule" being used over the extreme subintervals. 

For (x 2 - Xj) > .10, the maximum value of the function within the Interval 
Is found, and the approximation is performed over the coarse steps up to a point 
where the value of the function Is 20% of the maximum value. Beyond this point, 
the remaining coarse sublntexvals are each further divided Into 3 fine subintervals. 
The function Is evaluated at the endpoints of these fine subintervals. The Integral 
approximation Is completed using this finer step. 

For a given function f(x) evaluated at equally spaced (ax) points x Q , x^, 
x 2' x 3 # the qulntlc approximation for Integration Is given by: 

J 2 fW d. = (-f(x 0 ) + 13f(xj) + 

13f (x 2 ) - f(x 3 )) 

The "one-sided Simpson's Rule" is given by: 
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5,7,13 Subroutine QUITS 


Subroutine QUITS is entered when an error has been detected in the input 
or if an unreasonable quantity has been detected during execution. It prints out 
appropriate error statements and the contents of labeled common storage. 


5.7,14 Subroutine READSI 

Subroutine READSI is a modified version of subroutine READIN of Reference 
5. The name was changed to avoid a naming conflict with a subroutine by the 
same name in the ballistics module. 

This routine reads the input data for the TBL module in the $TBL namelist 
as described in Volume III, Section 2.10 o Subroutine READSI also checks the in- 
put data tor consistancy and prints out dianogstics. 

In addition to the above, this routine traps and/or reads any linkage data 
generated by the 0DE or TD2P modules. 
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5.7.15 Subroutine SEVAL 


Subroutine SEVAL, using certain defined Integral relationships, and a piece- 
wise cubic curve fit of specific heat (C p ) as a function of temperature (T), for a 
given Indicator value (Ind), evaluates the appropriate quantity: Pressure (P), en- 
tropy (S); enthalpy (H); temperature (T); specific heat (C p ) . 

For Ind = -1: P = P(T,S) 

0: S = S(T,P) 

1: H = H(T), C p = C p (t) 

2: T = T(H), C p = C p (t) 

3: T = T(S,P) 

The basic relations are: 

C p = A + BT + CP + DT 3 
H(T) = H. + f T C (t) dt 

1 md P 

T C nW P 

S(T) = S(T q ) + J T — E — dt - R In (y-) 

T o 
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5. 7 ,16 Subroutine START 

Subroutine START locates the sonic point (M=l) from the input Mach num- 
ber and x tables and then computes initial values for momentum and energy thick- 
nesses. 


The table of Mach numbers is first scanned to locate an exact value of 
M=l. If none is found, the sonic point location is determined from the following 
iteration procedure. 

1. Locate the i-th interval, in the table, which brackets M =1 and initially 
guess x g =x i# 

2. Using subroutine XNTERP, evaluate M g at x g . 

3. An absolute error comparison is made. If 

[ M g - l.| < TOLZME 

convergence is considered to be attained and appropriate flow quantities 
are evaluated at x . If unconverged, a secant method is used to cal- 
culate a new guess for x and steps 2 and 3 are repeated up to a max- 
imum of 50 iterations. g 

Using the flow quantities just computed, subroutine INTZET is used to 
evaluate integrals I ^ and I 2 for £= 1. The shape factor 5*/0 is then evaluated 
from equation (125 of Ref, 5). An iteration for the skin friction coefficient C f 
is entered (similar to that in subroutine BARPRO), and a value of momentum thick- 
ness q is computed. This value is assumed to be the initial value for 0 and </, 

at x, . 

3 


5 - 2.10 



5.7.17 Subroutine TBL 

Subroutine TBL initiates the execution of the TBL module by calling sub- 
routine DIRECT. 


5.7.18 Subroutine XNTERP 

Subroutine XNTERP, given tables of n points of independent variable (x^) 
and dependent variable (y^, returns point value (y) and first derivative (y 1 ) for a 
given value of x. 


A local quintic spline interpolation is defined as follows: The independent 
variable table is searched until the points bracketing x are found (x^ < x < * 1+1 ). 
Using the four surrounding points (x^, x^ x [+1 , x l+2 ) # coefficients of the 
parabola which pastes through the three leftmost points and that which passes 
through the three rightmost points are obtained. For the given value x, 

y L = c 1 + c 2 (x - x^j) + c 3 (x - x^) 2 

y R = c 4 +c 5 (x - x^) +c 6 (x - xf 

? L = c 2 + 2c 3 (x - 

? R = c 5 + 2c 6 Oc - x t ) 

A cubic weighting function ^and its derivative a' are defined as follows: 



U(x) 


x - x. 
x i+l “ x i 


ot (x) = 3U 2 - 2 IP 
' (x)= (6U - 6U 2 ) U' (x) 


Vhe dependent variable y and its derivative y' can now be evaluated. 
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5.7.19 Subroutine ZETAIT 


Subroutine ZETAIT calculates the shape parameter £ and boundary layer 
thicknesses a, fi and 5 * at a point x, for given values of 0 and 0 . 


For £ > 1 


£ = ( i //n ~ I 2'“ 1 3 
0 1 I, 


■) 


where 


For £<1 


6 


± 1 

e I^7T 

6 *_ / l/n " I 6~ I 7 x 

fi " ( "TTr ,r > 

0 a 4 a 5 


a = 


c = 


1 

n+l' 


6 = 




where 


C = 


± !& 

e i- 


i 

n+l 


An iteration procedure is used to evaluate £. 

1. An initial guess r is made 

Q 

2. For the value of £ ^1 or < 1) appropriate functions of p/p have been 

defined, and thr proper integrals are evaluated using subroutine INTZET 

3. £ is calculated by the appropriate formula a&ove 

4. A relative error comparison is made. If 


£ - £ 

| — t—2-I < TOLZET 

C g 

convergence is considered to be attained, if not converged, a form 
of Wegstein's method is useo to calculate a new guess for £ and 
steps 2 to 4 are repeated up to a maximum of 50 Iterations. g 

The boundary layer thicknesses 5 * and A are then calculated. 
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Appendix A - Conversions Aides 


The SPP code was developed under the rules and limitations of the SCOPE 
3.3 operating systems of CDC 6000 computers. The following notes are intended 
to aid the conversion of the SPP code to other computers and/or operating systems. 

1) Subroutine FIND of Section 5.1.6 is not the same routine as described 
in Section 5.6.5. On operating systems which poses a problem. 

It Is suggested that the name of the routine of 5.6.5 be changed to 

FINE. This routine is referred to in subroutines AC0MP, 0NED, PARTTL, 
and TRA CE. 

2) Subroutine ERROR of Section 5.4.8 is not the same as the routine de- 
scribed in 5.6.37. On operating systems where this poses a problem 
It is suggested that the routine of 5.4.8 should be changed to the name 
(external reference) of ERR0S. This routine (5.4.8) is only referred to 
in subroutine RE ADI N. 

3) A possible conflict in labeled common areas occurs in LINK41 and LINK50 
routines. Labeled common blocks by the name NAMEW occur in both 

of these overlays, however, these are not intended to be the same. It 
is suggested that the common block in subroutine 0DKENP be changed 
to the name NAMW to avoid this possible conflict. 

4) Contrary to the rules of F0RTRAN IV, EQUIVALENCE statements can 
sometime increase the size of labeled common areas. This can cause 
problems on some operating systems. In particular the labeled common 
block FNALBT in subroutine DERIV can be made to be 4 locations (cells) 
longer than defined in subroutine MAIN1D. If this causes problems, 
the length of this common block should be extended by 4 dummy loca- 
tions (which are not used) in subroutines MAIN ID, FLU, and IAUX. 

Other minor corrections to the SPP code which will update the SPP code to 
the 1.2 version are: 

a) Moving the labeled common block D0L00PS into the 0DK routine. 

b) Adding the variable CPS0L to the $TD2 namelist. 

c) Adding the test at approximately line 65 of subroutine TD2P of 

IF (FTD2P .EQ. 0 . 0 ) CPS0L = 0.0 

d) Changing the test is subroutine AGP from 

IF (FTD2P .EQ.0) G0 T0 12 
to: IF(CPS0L .EQ.0) G0 T0 12 
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The conversion of the SPP code to the Univac 1108, EXEC 8 operation system 
took approximately one man day of effort. However, it is recognized by long and 
painful years of experience that the authors of a particular code have much less 
trouble converting codes than those unfamiliar with the particular program in 
question. Therefore, questions concerning the SPP code should be directed to the 
Air Force contracting officier in charge of the SPP code and responses will be made 
through him. 

Other notes concerning conversions of the SPP code are as follows for 
Univac and IBM users. 

I) The cards with 0VERLAY (SPP,i,j) should be removed. 

II) The cards with PR0GRAM LINKij should be changed to SUBR0UHNE 
LINKiJ and a RETURN card should be added just before the END card. 

III) The calls to 0VERLAY (3LSPP,l,j) should be changed to CALL LINKij. 

It is hoped that these brief notes will substantially aid in the conversion of 
the SPP code to other computer and/or operating systems for which the SPP code 
was not designed for. 
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